home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / gnu / gcc_260.zip / gcc_260 / real.c < prev    next >
C/C++ Source or Header  |  1994-06-28  |  123KB  |  5,970 lines
  1. /* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
  2.    and support for XFmode IEEE extended real floating point arithmetic.
  3.    Copyright (C) 1993, 1994 Free Software Foundation, Inc.
  4.    Contributed by Stephen L. Moshier (moshier@world.std.com).
  5.  
  6. This file is part of GNU CC.
  7.  
  8. GNU CC is free software; you can redistribute it and/or modify
  9. it under the terms of the GNU General Public License as published by
  10. the Free Software Foundation; either version 2, or (at your option)
  11. any later version.
  12.  
  13. GNU CC is distributed in the hope that it will be useful,
  14. but WITHOUT ANY WARRANTY; without even the implied warranty of
  15. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16. GNU General Public License for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with GNU CC; see the file COPYING.  If not, write to
  20. the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
  21.  
  22. #include <stdio.h>
  23. #include <errno.h>
  24. #include "config.h"
  25. #include "tree.h"
  26.  
  27. #ifndef errno
  28. extern int errno;
  29. #endif
  30.  
  31. /* To enable support of XFmode extended real floating point, define
  32. LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
  33.  
  34. To support cross compilation between IEEE, VAX and IBM floating
  35. point formats, define REAL_ARITHMETIC in the tm.h file.
  36.  
  37. In either case the machine files (tm.h) must not contain any code
  38. that tries to use host floating point arithmetic to convert
  39. REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
  40. etc.  In cross-compile situations a REAL_VALUE_TYPE may not
  41. be intelligible to the host computer's native arithmetic.
  42.  
  43. The emulator defaults to the host's floating point format so that
  44. its decimal conversion functions can be used if desired (see
  45. real.h).
  46.  
  47. The first part of this file interfaces gcc to ieee.c, which is a
  48. floating point arithmetic suite that was not written with gcc in
  49. mind.  The interface is followed by ieee.c itself and related
  50. items. Avoid changing ieee.c unless you have suitable test
  51. programs available.  A special version of the PARANOIA floating
  52. point arithmetic tester, modified for this purpose, can be found
  53. on usc.edu : /pub/C-numanal/ieeetest.zoo.  Some tutorial
  54. information on ieee.c is given in my book: S. L. Moshier,
  55. _Methods and Programs for Mathematical Functions_, Prentice-Hall
  56. or Simon & Schuster Int'l, 1989.  A library of XFmode elementary
  57. transcendental functions can be obtained by ftp from
  58. research.att.com: netlib/cephes/ldouble.shar.Z  */
  59.  
  60. /* Type of computer arithmetic.
  61.    Only one of DEC, IBM, MIEEE, IBMPC, or UNK should get defined.
  62.  
  63.    `MIEEE' refers generically to big-endian IEEE floating-point data
  64.    structure.  This definition should work in SFmode `float' type and
  65.    DFmode `double' type on virtually all big-endian IEEE machines.
  66.    If LONG_DOUBLE_TYPE_SIZE has been defined to be 96, then MIEEE
  67.    also invokes the particular XFmode (`long double' type) data
  68.    structure used by the Motorola 680x0 series processors.
  69.  
  70.    `IBMPC' refers generally to little-endian IEEE machines. In this
  71.    case, if LONG_DOUBLE_TYPE_SIZE has been defined to be 96, then
  72.    IBMPC also invokes the particular XFmode `long double' data
  73.    structure used by the Intel 80x86 series processors.
  74.  
  75.    `DEC' refers specifically to the Digital Equipment Corp PDP-11
  76.    and VAX floating point data structure.  This model currently
  77.    supports no type wider than DFmode.
  78.  
  79.    `IBM' refers specifically to the IBM System/370 and compatible
  80.    floating point data structure.  This model currently supports
  81.    no type wider than DFmode.  The IBM conversions were contributed by
  82.    frank@atom.ansto.gov.au (Frank Crawford).
  83.  
  84.    If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
  85.    then `long double' and `double' are both implemented, but they
  86.    both mean DFmode.  In this case, the software floating-point
  87.    support available here is activated by writing
  88.       #define REAL_ARITHMETIC
  89.    in tm.h. 
  90.  
  91.    The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
  92.    and may deactivate XFmode since `long double' is used to refer
  93.    to both modes.
  94.  
  95.    The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
  96.    contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
  97.    separate the floating point unit's endian-ness from that of
  98.    the integer addressing.  This permits one to define a big-endian
  99.    FPU on a little-endian machine (e.g., ARM).  An extension to
  100.    BYTES_BIG_ENDIAN may be required for some machines in the future.
  101.    These optional macros may be defined in tm.h.  In real.h, they
  102.    default to WORDS_BIG_ENDIAN, etc., so there is no need to define
  103.    them for any normal host or target machine on which the floats
  104.    and the integers have the same endian-ness.   */
  105.  
  106.  
  107. /* The following converts gcc macros into the ones used by this file.  */
  108.  
  109. /* REAL_ARITHMETIC defined means that macros in real.h are
  110.    defined to call emulator functions.  */
  111. #ifdef REAL_ARITHMETIC
  112.  
  113. #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
  114. /* PDP-11, Pro350, VAX: */
  115. #define DEC 1
  116. #else /* it's not VAX */
  117. #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
  118. /* IBM System/370 style */
  119. #define IBM 1
  120. #else /* it's also not an IBM */
  121. #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
  122. #if FLOAT_WORDS_BIG_ENDIAN
  123. /* Motorola IEEE, high order words come first (Sun workstation): */
  124. #define MIEEE 1
  125. #else /* not big-endian */
  126. /* Intel IEEE, low order words come first:
  127.  */
  128. #define IBMPC 1
  129. #endif /*  big-endian */
  130. #else /* it's not IEEE either */
  131. /* UNKnown arithmetic.  We don't support this and can't go on. */
  132. unknown arithmetic type
  133. #define UNK 1
  134. #endif /* not IEEE */
  135. #endif /* not IBM */
  136. #endif /* not VAX */
  137.  
  138. #else
  139. /* REAL_ARITHMETIC not defined means that the *host's* data
  140.    structure will be used.  It may differ by endian-ness from the
  141.    target machine's structure and will get its ends swapped
  142.    accordingly (but not here).  Probably only the decimal <-> binary
  143.    functions in this file will actually be used in this case.  */
  144.  
  145. #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
  146. #define DEC 1
  147. #else /* it's not VAX */
  148. #if HOST_FLOAT_FORMAT == IBM_FLOAT_FORMAT
  149. /* IBM System/370 style */
  150. #define IBM 1
  151. #else /* it's also not an IBM */
  152. #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
  153. #if HOST_FLOAT_WORDS_BIG_ENDIAN
  154. #define MIEEE 1
  155. #else /* not big-endian */
  156. #define IBMPC 1
  157. #endif /*  big-endian */
  158. #else /* it's not IEEE either */
  159. unknown arithmetic type
  160. #define UNK 1
  161. #endif /* not IEEE */
  162. #endif /* not IBM */
  163. #endif /* not VAX */
  164.  
  165. #endif /* REAL_ARITHMETIC not defined */
  166.  
  167. /* Define INFINITY for support of infinity.
  168.    Define NANS for support of Not-a-Number's (NaN's).  */
  169. #if !defined(DEC) && !defined(IBM)
  170. #define INFINITY
  171. #define NANS
  172. #endif
  173.  
  174. /* Support of NaNs requires support of infinity. */
  175. #ifdef NANS
  176. #ifndef INFINITY
  177. #define INFINITY
  178. #endif
  179. #endif
  180.  
  181. /* Find a host integer type that is at least 16 bits wide,
  182.    and another type at least twice whatever that size is. */
  183.  
  184. #if HOST_BITS_PER_CHAR >= 16
  185. #define EMUSHORT char
  186. #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
  187. #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
  188. #else
  189. #if HOST_BITS_PER_SHORT >= 16
  190. #define EMUSHORT short
  191. #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
  192. #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
  193. #else
  194. #if HOST_BITS_PER_INT >= 16
  195. #define EMUSHORT int
  196. #define EMUSHORT_SIZE HOST_BITS_PER_INT
  197. #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
  198. #else
  199. #if HOST_BITS_PER_LONG >= 16
  200. #define EMUSHORT long
  201. #define EMUSHORT_SIZE HOST_BITS_PER_LONG
  202. #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
  203. #else
  204. /*  You will have to modify this program to have a smaller unit size. */
  205. #define EMU_NON_COMPILE
  206. #endif
  207. #endif
  208. #endif
  209. #endif
  210.  
  211. #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
  212. #define EMULONG short
  213. #else
  214. #if HOST_BITS_PER_INT >= EMULONG_SIZE
  215. #define EMULONG int
  216. #else
  217. #if HOST_BITS_PER_LONG >= EMULONG_SIZE
  218. #define EMULONG long
  219. #else
  220. #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
  221. #define EMULONG long long int
  222. #else
  223. /*  You will have to modify this program to have a smaller unit size. */
  224. #define EMU_NON_COMPILE
  225. #endif
  226. #endif
  227. #endif
  228. #endif
  229.  
  230.  
  231. /* The host interface doesn't work if no 16-bit size exists. */
  232. #if EMUSHORT_SIZE != 16
  233. #define EMU_NON_COMPILE
  234. #endif
  235.  
  236. /* OK to continue compilation. */
  237. #ifndef EMU_NON_COMPILE
  238.  
  239. /* Construct macros to translate between REAL_VALUE_TYPE and e type.
  240.    In GET_REAL and PUT_REAL, r and e are pointers.
  241.    A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
  242.    in memory, with no holes.  */
  243.  
  244. #if LONG_DOUBLE_TYPE_SIZE == 96
  245. /* Number of 16 bit words in external e type format */
  246. #define NE 6
  247. #define MAXDECEXP 4932
  248. #define MINDECEXP -4956
  249. #define GET_REAL(r,e) bcopy (r, e, 2*NE)
  250. #define PUT_REAL(e,r) bcopy (e, r, 2*NE)
  251. #else /* no XFmode */
  252. #if LONG_DOUBLE_TYPE_SIZE == 128
  253. #define NE 10
  254. #define MAXDECEXP 4932
  255. #define MINDECEXP -4977
  256. #define GET_REAL(r,e) bcopy (r, e, 2*NE)
  257. #define PUT_REAL(e,r) bcopy (e, r, 2*NE)
  258. #else
  259. #define NE 6
  260. #define MAXDECEXP 4932
  261. #define MINDECEXP -4956
  262. #ifdef REAL_ARITHMETIC
  263. /* Emulator uses target format internally
  264.    but host stores it in host endian-ness. */
  265.  
  266. #if HOST_FLOAT_WORDS_BIG_ENDIAN == FLOAT_WORDS_BIG_ENDIAN
  267. #define GET_REAL(r,e) e53toe ((unsigned EMUSHORT*) (r), (e))
  268. #define PUT_REAL(e,r) etoe53 ((e), (unsigned EMUSHORT *) (r))
  269.  
  270. #else /* endian-ness differs */
  271. /* emulator uses target endian-ness internally */
  272. #define GET_REAL(r,e)        \
  273. do { unsigned EMUSHORT w[4];    \
  274.  w[3] = ((EMUSHORT *) r)[0];    \
  275.  w[2] = ((EMUSHORT *) r)[1];    \
  276.  w[1] = ((EMUSHORT *) r)[2];    \
  277.  w[0] = ((EMUSHORT *) r)[3];    \
  278.  e53toe (w, (e)); } while (0)
  279.  
  280. #define PUT_REAL(e,r)        \
  281. do { unsigned EMUSHORT w[4];    \
  282.  etoe53 ((e), w);        \
  283.  *((EMUSHORT *) r) = w[3];    \
  284.  *((EMUSHORT *) r + 1) = w[2];    \
  285.  *((EMUSHORT *) r + 2) = w[1];    \
  286.  *((EMUSHORT *) r + 3) = w[0]; } while (0)
  287.  
  288. #endif /* endian-ness differs */
  289.  
  290. #else /* not REAL_ARITHMETIC */
  291.  
  292. /* emulator uses host format */
  293. #define GET_REAL(r,e) e53toe ((unsigned EMUSHORT *) (r), (e))
  294. #define PUT_REAL(e,r) etoe53 ((e), (unsigned EMUSHORT *) (r))
  295.  
  296. #endif /* not REAL_ARITHMETIC */
  297. #endif /* not TFmode */
  298. #endif /* no XFmode */
  299.  
  300.  
  301. /* Number of 16 bit words in internal format */
  302. #define NI (NE+3)
  303.  
  304. /* Array offset to exponent */
  305. #define E 1
  306.  
  307. /* Array offset to high guard word */
  308. #define M 2
  309.  
  310. /* Number of bits of precision */
  311. #define NBITS ((NI-4)*16)
  312.  
  313. /* Maximum number of decimal digits in ASCII conversion
  314.  * = NBITS*log10(2)
  315.  */
  316. #define NDEC (NBITS*8/27)
  317.  
  318. /* The exponent of 1.0 */
  319. #define EXONE (0x3fff)
  320.  
  321. extern int extra_warnings;
  322. extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
  323. extern unsigned EMUSHORT elog2[], esqrt2[];
  324.  
  325. static void endian    PROTO((unsigned EMUSHORT *, long *,
  326.                    enum machine_mode));
  327. static void eclear    PROTO((unsigned EMUSHORT *));
  328. static void emov    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  329. static void eabs    PROTO((unsigned EMUSHORT *));
  330. static void eneg    PROTO((unsigned EMUSHORT *));
  331. static int eisneg    PROTO((unsigned EMUSHORT *));
  332. static int eisinf    PROTO((unsigned EMUSHORT *));
  333. static int eisnan    PROTO((unsigned EMUSHORT *));
  334. static void einfin    PROTO((unsigned EMUSHORT *));
  335. static void enan    PROTO((unsigned EMUSHORT *, int));
  336. static void emovi    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  337. static void emovo    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  338. static void ecleaz    PROTO((unsigned EMUSHORT *));
  339. static void ecleazs    PROTO((unsigned EMUSHORT *));
  340. static void emovz    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  341. static void einan    PROTO((unsigned EMUSHORT *));
  342. static int eiisnan    PROTO((unsigned EMUSHORT *));
  343. static int eiisneg    PROTO((unsigned EMUSHORT *));
  344. static void eiinfin    PROTO((unsigned EMUSHORT *));
  345. static int eiisinf    PROTO((unsigned EMUSHORT *));
  346. static int ecmpm    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  347. static void eshdn1    PROTO((unsigned EMUSHORT *));
  348. static void eshup1    PROTO((unsigned EMUSHORT *));
  349. static void eshdn8    PROTO((unsigned EMUSHORT *));
  350. static void eshup8    PROTO((unsigned EMUSHORT *));
  351. static void eshup6    PROTO((unsigned EMUSHORT *));
  352. static void eshdn6    PROTO((unsigned EMUSHORT *));
  353. static void eaddm    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  354. static void esubm    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  355. static void m16m    PROTO((unsigned int, unsigned short *,
  356.                    unsigned short *));
  357. static int edivm    PROTO((unsigned short *, unsigned short *));
  358. static int emulm    PROTO((unsigned short *, unsigned short *));
  359. static void emdnorm    PROTO((unsigned EMUSHORT *, int, int, EMULONG, int));
  360. static void esub    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  361.                    unsigned EMUSHORT *));
  362. static void eadd    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  363.                    unsigned EMUSHORT *));
  364. static void eadd1    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  365.                    unsigned EMUSHORT *));
  366. static void ediv    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  367.                    unsigned EMUSHORT *));
  368. static void emul    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  369.                    unsigned EMUSHORT *));
  370. static void e53toe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  371. static void e64toe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  372. static void e113toe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  373. static void e24toe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  374. static void etoe113    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  375. static void toe113    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  376. static void etoe64    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  377. static void toe64    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  378. static void etoe53    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  379. static void toe53    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  380. static void etoe24    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  381. static void toe24    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  382. static int ecmp        PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  383. static void eround    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  384. static void ltoe    PROTO((HOST_WIDE_INT *, unsigned EMUSHORT *));
  385. static void ultoe    PROTO((unsigned HOST_WIDE_INT *, unsigned EMUSHORT *));
  386. static void eifrac    PROTO((unsigned EMUSHORT *, HOST_WIDE_INT *,
  387.                    unsigned EMUSHORT *));
  388. static void euifrac    PROTO((unsigned EMUSHORT *, unsigned HOST_WIDE_INT *,
  389.                    unsigned EMUSHORT *));
  390. static int eshift    PROTO((unsigned EMUSHORT *, int));
  391. static int enormlz    PROTO((unsigned EMUSHORT *));
  392. static void e24toasc    PROTO((unsigned EMUSHORT *, char *, int));
  393. static void e53toasc    PROTO((unsigned EMUSHORT *, char *, int));
  394. static void e64toasc    PROTO((unsigned EMUSHORT *, char *, int));
  395. static void e113toasc    PROTO((unsigned EMUSHORT *, char *, int));
  396. static void etoasc    PROTO((unsigned EMUSHORT *, char *, int));
  397. static void asctoe24    PROTO((char *, unsigned EMUSHORT *));
  398. static void asctoe53    PROTO((char *, unsigned EMUSHORT *));
  399. static void asctoe64    PROTO((char *, unsigned EMUSHORT *));
  400. static void asctoe113    PROTO((char *, unsigned EMUSHORT *));
  401. static void asctoe    PROTO((char *, unsigned EMUSHORT *));
  402. static void asctoeg    PROTO((char *, unsigned EMUSHORT *, int));
  403. static void efloor    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  404. static void efrexp    PROTO((unsigned EMUSHORT *, int *,
  405.                    unsigned EMUSHORT *));
  406. static void eldexp    PROTO((unsigned EMUSHORT *, int, unsigned EMUSHORT *));
  407. static void eremain    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  408.                    unsigned EMUSHORT *));
  409. static void eiremain    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  410. static void mtherr    PROTO((char *, int));
  411. static void dectoe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  412. static void etodec    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  413. static void todec    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  414. static void ibmtoe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  415.                    enum machine_mode));
  416. static void etoibm    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  417.                    enum machine_mode));
  418. static void toibm    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
  419.                    enum machine_mode));
  420. static void make_nan    PROTO((unsigned EMUSHORT *, int, enum machine_mode));
  421. static void uditoe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  422. static void ditoe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  423. static void etoudi    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  424. static void etodi    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  425. static void esqrt    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
  426.  
  427. /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
  428.    swapping ends if required, into output array of longs.  The
  429.    result is normally passed to fprintf by the ASM_OUTPUT_ macros.   */
  430.  
  431. static void 
  432. endian (e, x, mode)
  433.      unsigned EMUSHORT e[];
  434.      long x[];
  435.      enum machine_mode mode;
  436. {
  437.   unsigned long th, t;
  438.  
  439. #if FLOAT_WORDS_BIG_ENDIAN
  440.   switch (mode)
  441.     {
  442.  
  443.     case TFmode:
  444.       /* Swap halfwords in the fourth long. */
  445.       th = (unsigned long) e[6] & 0xffff;
  446.       t = (unsigned long) e[7] & 0xffff;
  447.       t |= th << 16;
  448.       x[3] = (long) t;
  449.  
  450.     case XFmode:
  451.  
  452.       /* Swap halfwords in the third long. */
  453.       th = (unsigned long) e[4] & 0xffff;
  454.       t = (unsigned long) e[5] & 0xffff;
  455.       t |= th << 16;
  456.       x[2] = (long) t;
  457.       /* fall into the double case */
  458.  
  459.     case DFmode:
  460.  
  461.       /* swap halfwords in the second word */
  462.       th = (unsigned long) e[2] & 0xffff;
  463.       t = (unsigned long) e[3] & 0xffff;
  464.       t |= th << 16;
  465.       x[1] = (long) t;
  466.       /* fall into the float case */
  467.  
  468.     case HFmode:
  469.     case SFmode:
  470.  
  471.       /* swap halfwords in the first word */
  472.       th = (unsigned long) e[0] & 0xffff;
  473.       t = (unsigned long) e[1] & 0xffff;
  474.       t |= th << 16;
  475.       x[0] = t;
  476.       break;
  477.  
  478.     default:
  479.       abort ();
  480.     }
  481.  
  482. #else
  483.  
  484.   /* Pack the output array without swapping. */
  485.  
  486.   switch (mode)
  487.     {
  488.  
  489.     case TFmode:
  490.  
  491.       /* Pack the fourth long. */
  492.       th = (unsigned long) e[7] & 0xffff;
  493.       t = (unsigned long) e[6] & 0xffff;
  494.       t |= th << 16;
  495.       x[3] = (long) t;
  496.  
  497.     case XFmode:
  498.  
  499.       /* Pack the third long.
  500.      Each element of the input REAL_VALUE_TYPE array has 16 useful bits
  501.      in it.  */
  502.       th = (unsigned long) e[5] & 0xffff;
  503.       t = (unsigned long) e[4] & 0xffff;
  504.       t |= th << 16;
  505.       x[2] = (long) t;
  506.       /* fall into the double case */
  507.  
  508.     case DFmode:
  509.  
  510.       /* pack the second long */
  511.       th = (unsigned long) e[3] & 0xffff;
  512.       t = (unsigned long) e[2] & 0xffff;
  513.       t |= th << 16;
  514.       x[1] = (long) t;
  515.       /* fall into the float case */
  516.  
  517.     case HFmode:
  518.     case SFmode:
  519.  
  520.       /* pack the first long */
  521.       th = (unsigned long) e[1] & 0xffff;
  522.       t = (unsigned long) e[0] & 0xffff;
  523.       t |= th << 16;
  524.       x[0] = t;
  525.       break;
  526.  
  527.     default:
  528.       abort ();
  529.     }
  530.  
  531. #endif
  532. }
  533.  
  534.  
  535. /* This is the implementation of the REAL_ARITHMETIC macro.  */
  536.  
  537. void 
  538. earith (value, icode, r1, r2)
  539.      REAL_VALUE_TYPE *value;
  540.      int icode;
  541.      REAL_VALUE_TYPE *r1;
  542.      REAL_VALUE_TYPE *r2;
  543. {
  544.   unsigned EMUSHORT d1[NE], d2[NE], v[NE];
  545.   enum tree_code code;
  546.  
  547.   GET_REAL (r1, d1);
  548.   GET_REAL (r2, d2);
  549. #ifdef NANS
  550. /*  Return NaN input back to the caller. */
  551.   if (eisnan (d1))
  552.     {
  553.       PUT_REAL (d1, value);
  554.       return;
  555.     }
  556.   if (eisnan (d2))
  557.     {
  558.       PUT_REAL (d2, value);
  559.       return;
  560.     }
  561. #endif
  562.   code = (enum tree_code) icode;
  563.   switch (code)
  564.     {
  565.     case PLUS_EXPR:
  566.       eadd (d2, d1, v);
  567.       break;
  568.  
  569.     case MINUS_EXPR:
  570.       esub (d2, d1, v);        /* d1 - d2 */
  571.       break;
  572.  
  573.     case MULT_EXPR:
  574.       emul (d2, d1, v);
  575.       break;
  576.  
  577.     case RDIV_EXPR:
  578. #ifndef REAL_INFINITY
  579.       if (ecmp (d2, ezero) == 0)
  580.     {
  581. #ifdef NANS
  582.     enan (v, eisneg (d1) ^ eisneg (d2));
  583.     break;
  584. #else
  585.     abort ();
  586. #endif
  587.     }
  588. #endif
  589.       ediv (d2, d1, v);    /* d1/d2 */
  590.       break;
  591.  
  592.     case MIN_EXPR:        /* min (d1,d2) */
  593.       if (ecmp (d1, d2) < 0)
  594.     emov (d1, v);
  595.       else
  596.     emov (d2, v);
  597.       break;
  598.  
  599.     case MAX_EXPR:        /* max (d1,d2) */
  600.       if (ecmp (d1, d2) > 0)
  601.     emov (d1, v);
  602.       else
  603.     emov (d2, v);
  604.       break;
  605.     default:
  606.       emov (ezero, v);
  607.       break;
  608.     }
  609. PUT_REAL (v, value);
  610. }
  611.  
  612.  
  613. /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
  614.    implements REAL_VALUE_RNDZINT (x) (etrunci (x)).  */
  615.  
  616. REAL_VALUE_TYPE 
  617. etrunci (x)
  618.      REAL_VALUE_TYPE x;
  619. {
  620.   unsigned EMUSHORT f[NE], g[NE];
  621.   REAL_VALUE_TYPE r;
  622.   HOST_WIDE_INT l;
  623.  
  624.   GET_REAL (&x, g);
  625. #ifdef NANS
  626.   if (eisnan (g))
  627.     return (x);
  628. #endif
  629.   eifrac (g, &l, f);
  630.   ltoe (&l, g);
  631.   PUT_REAL (g, &r);
  632.   return (r);
  633. }
  634.  
  635.  
  636. /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
  637.    implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)).  */
  638.  
  639. REAL_VALUE_TYPE 
  640. etruncui (x)
  641.      REAL_VALUE_TYPE x;
  642. {
  643.   unsigned EMUSHORT f[NE], g[NE];
  644.   REAL_VALUE_TYPE r;
  645.   unsigned HOST_WIDE_INT l;
  646.  
  647.   GET_REAL (&x, g);
  648. #ifdef NANS
  649.   if (eisnan (g))
  650.     return (x);
  651. #endif
  652.   euifrac (g, &l, f);
  653.   ultoe (&l, g);
  654.   PUT_REAL (g, &r);
  655.   return (r);
  656. }
  657.  
  658.  
  659. /* This is the REAL_VALUE_ATOF function.  It converts a decimal string to
  660.    binary, rounding off as indicated by the machine_mode argument.  Then it
  661.    promotes the rounded value to REAL_VALUE_TYPE.  */
  662.  
  663. REAL_VALUE_TYPE 
  664. ereal_atof (s, t)
  665.      char *s;
  666.      enum machine_mode t;
  667. {
  668.   unsigned EMUSHORT tem[NE], e[NE];
  669.   REAL_VALUE_TYPE r;
  670.  
  671.   switch (t)
  672.     {
  673.     case HFmode:
  674.     case SFmode:
  675.       asctoe24 (s, tem);
  676.       e24toe (tem, e);
  677.       break;
  678.     case DFmode:
  679.       asctoe53 (s, tem);
  680.       e53toe (tem, e);
  681.       break;
  682.     case XFmode:
  683.       asctoe64 (s, tem);
  684.       e64toe (tem, e);
  685.       break;
  686.     case TFmode:
  687.       asctoe113 (s, tem);
  688.       e113toe (tem, e);
  689.       break;
  690.     default:
  691.       asctoe (s, e);
  692.     }
  693.   PUT_REAL (e, &r);
  694.   return (r);
  695. }
  696.  
  697.  
  698. /* Expansion of REAL_NEGATE.  */
  699.  
  700. REAL_VALUE_TYPE 
  701. ereal_negate (x)
  702.      REAL_VALUE_TYPE x;
  703. {
  704.   unsigned EMUSHORT e[NE];
  705.   REAL_VALUE_TYPE r;
  706.  
  707.   GET_REAL (&x, e);
  708.   eneg (e);
  709.   PUT_REAL (e, &r);
  710.   return (r);
  711. }
  712.  
  713.  
  714. /* Round real toward zero to HOST_WIDE_INT;
  715.    implements REAL_VALUE_FIX (x).  */
  716.  
  717. HOST_WIDE_INT
  718. efixi (x)
  719.      REAL_VALUE_TYPE x;
  720. {
  721.   unsigned EMUSHORT f[NE], g[NE];
  722.   HOST_WIDE_INT l;
  723.  
  724.   GET_REAL (&x, f);
  725. #ifdef NANS
  726.   if (eisnan (f))
  727.     {
  728.       warning ("conversion from NaN to int");
  729.       return (-1);
  730.     }
  731. #endif
  732.   eifrac (f, &l, g);
  733.   return l;
  734. }
  735.  
  736. /* Round real toward zero to unsigned HOST_WIDE_INT
  737.    implements  REAL_VALUE_UNSIGNED_FIX (x).
  738.    Negative input returns zero.  */
  739.  
  740. unsigned HOST_WIDE_INT
  741. efixui (x)
  742.      REAL_VALUE_TYPE x;
  743. {
  744.   unsigned EMUSHORT f[NE], g[NE];
  745.   unsigned HOST_WIDE_INT l;
  746.  
  747.   GET_REAL (&x, f);
  748. #ifdef NANS
  749.   if (eisnan (f))
  750.     {
  751.       warning ("conversion from NaN to unsigned int");
  752.       return (-1);
  753.     }
  754. #endif
  755.   euifrac (f, &l, g);
  756.   return l;
  757. }
  758.  
  759.  
  760. /* REAL_VALUE_FROM_INT macro.  */
  761.  
  762. void 
  763. ereal_from_int (d, i, j)
  764.      REAL_VALUE_TYPE *d;
  765.      HOST_WIDE_INT i, j;
  766. {
  767.   unsigned EMUSHORT df[NE], dg[NE];
  768.   HOST_WIDE_INT low, high;
  769.   int sign;
  770.  
  771.   sign = 0;
  772.   low = i;
  773.   if ((high = j) < 0)
  774.     {
  775.       sign = 1;
  776.       /* complement and add 1 */
  777.       high = ~high;
  778.       if (low)
  779.     low = -low;
  780.       else
  781.     high += 1;
  782.     }
  783.   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  784.   ultoe ((unsigned HOST_WIDE_INT *) &high, dg);
  785.   emul (dg, df, dg);
  786.   ultoe ((unsigned HOST_WIDE_INT *) &low, df);
  787.   eadd (df, dg, dg);
  788.   if (sign)
  789.     eneg (dg);
  790.   PUT_REAL (dg, d);
  791. }
  792.  
  793.  
  794. /* REAL_VALUE_FROM_UNSIGNED_INT macro.   */
  795.  
  796. void 
  797. ereal_from_uint (d, i, j)
  798.      REAL_VALUE_TYPE *d;
  799.      unsigned HOST_WIDE_INT i, j;
  800. {
  801.   unsigned EMUSHORT df[NE], dg[NE];
  802.   unsigned HOST_WIDE_INT low, high;
  803.  
  804.   low = i;
  805.   high = j;
  806.   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  807.   ultoe (&high, dg);
  808.   emul (dg, df, dg);
  809.   ultoe (&low, df);
  810.   eadd (df, dg, dg);
  811.   PUT_REAL (dg, d);
  812. }
  813.  
  814.  
  815. /* REAL_VALUE_TO_INT macro.  */
  816.  
  817. void 
  818. ereal_to_int (low, high, rr)
  819.      HOST_WIDE_INT *low, *high;
  820.      REAL_VALUE_TYPE rr;
  821. {
  822.   unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
  823.   int s;
  824.  
  825.   GET_REAL (&rr, d);
  826. #ifdef NANS
  827.   if (eisnan (d))
  828.     {
  829.       warning ("conversion from NaN to int");
  830.       *low = -1;
  831.       *high = -1;
  832.       return;
  833.     }
  834. #endif
  835.   /* convert positive value */
  836.   s = 0;
  837.   if (eisneg (d))
  838.     {
  839.       eneg (d);
  840.       s = 1;
  841.     }
  842.   eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  843.   ediv (df, d, dg);        /* dg = d / 2^32 is the high word */
  844.   euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);
  845.   emul (df, dh, dg);        /* fractional part is the low word */
  846.   euifrac (dg, (unsigned HOST_WIDE_INT *)low, dh);
  847.   if (s)
  848.     {
  849.       /* complement and add 1 */
  850.       *high = ~(*high);
  851.       if (*low)
  852.     *low = -(*low);
  853.       else
  854.     *high += 1;
  855.     }
  856. }
  857.  
  858.  
  859. /* REAL_VALUE_LDEXP macro.  */
  860.  
  861. REAL_VALUE_TYPE
  862. ereal_ldexp (x, n)
  863.      REAL_VALUE_TYPE x;
  864.      int n;
  865. {
  866.   unsigned EMUSHORT e[NE], y[NE];
  867.   REAL_VALUE_TYPE r;
  868.  
  869.   GET_REAL (&x, e);
  870. #ifdef NANS
  871.   if (eisnan (e))
  872.     return (x);
  873. #endif
  874.   eldexp (e, n, y);
  875.   PUT_REAL (y, &r);
  876.   return (r);
  877. }
  878.  
  879. /* These routines are conditionally compiled because functions
  880.    of the same names may be defined in fold-const.c.  */
  881.  
  882. #ifdef REAL_ARITHMETIC
  883.  
  884. /* Check for infinity in a REAL_VALUE_TYPE. */
  885.  
  886. int
  887. target_isinf (x)
  888.      REAL_VALUE_TYPE x;
  889. {
  890.   unsigned EMUSHORT e[NE];
  891.  
  892. #ifdef INFINITY
  893.   GET_REAL (&x, e);
  894.   return (eisinf (e));
  895. #else
  896.   return 0;
  897. #endif
  898. }
  899.  
  900.  
  901. /* Check whether a REAL_VALUE_TYPE item is a NaN. */
  902.  
  903. int
  904. target_isnan (x)
  905.      REAL_VALUE_TYPE x;
  906. {
  907.   unsigned EMUSHORT e[NE];
  908.  
  909. #ifdef NANS
  910.   GET_REAL (&x, e);
  911.   return (eisnan (e));
  912. #else
  913.   return (0);
  914. #endif
  915. }
  916.  
  917.  
  918. /* Check for a negative REAL_VALUE_TYPE number.
  919.    This just checks the sign bit, so that -0 counts as negative. */
  920.  
  921. int
  922. target_negative (x)
  923.      REAL_VALUE_TYPE x;
  924. {
  925.   return ereal_isneg (x);
  926. }
  927.  
  928. /* Expansion of REAL_VALUE_TRUNCATE.
  929.    The result is in floating point, rounded to nearest or even.  */
  930.  
  931. REAL_VALUE_TYPE
  932. real_value_truncate (mode, arg)
  933.      enum machine_mode mode;
  934.      REAL_VALUE_TYPE arg;
  935. {
  936.   unsigned EMUSHORT e[NE], t[NE];
  937.   REAL_VALUE_TYPE r;
  938.  
  939.   GET_REAL (&arg, e);
  940. #ifdef NANS
  941.   if (eisnan (e))
  942.     return (arg);
  943. #endif
  944.   eclear (t);
  945.   switch (mode)
  946.     {
  947.     case TFmode:
  948.       etoe113 (e, t);
  949.       e113toe (t, t);
  950.       break;
  951.  
  952.     case XFmode:
  953.       etoe64 (e, t);
  954.       e64toe (t, t);
  955.       break;
  956.  
  957.     case DFmode:
  958.       etoe53 (e, t);
  959.       e53toe (t, t);
  960.       break;
  961.  
  962.     case HFmode:
  963.     case SFmode:
  964.       etoe24 (e, t);
  965.       e24toe (t, t);
  966.       break;
  967.  
  968.     case SImode:
  969.       r = etrunci (arg);
  970.       return (r);
  971.  
  972.     /* If an unsupported type was requested, presume that
  973.        the machine files know something useful to do with
  974.        the unmodified value.  */
  975.  
  976.     default:
  977.       return (arg);
  978.     }
  979.   PUT_REAL (t, &r);
  980.   return (r);
  981. }
  982.  
  983. #endif /* REAL_ARITHMETIC defined */
  984.  
  985. /* Used for debugging--print the value of R in human-readable format
  986.    on stderr.  */
  987.  
  988. void
  989. debug_real (r)
  990.      REAL_VALUE_TYPE r;
  991. {
  992.   char dstr[30];
  993.  
  994.   REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
  995.   fprintf (stderr, "%s", dstr);
  996. }  
  997.  
  998.  
  999. /* Target values are arrays of host longs. A long is guaranteed
  1000.    to be at least 32 bits wide. */
  1001.  
  1002. /* 128-bit long double */
  1003.  
  1004. void 
  1005. etartdouble (r, l)
  1006.      REAL_VALUE_TYPE r;
  1007.      long l[];
  1008. {
  1009.   unsigned EMUSHORT e[NE];
  1010.  
  1011.   GET_REAL (&r, e);
  1012.   etoe113 (e, e);
  1013.   endian (e, l, TFmode);
  1014. }
  1015.  
  1016. /* 80-bit long double */
  1017.  
  1018. void 
  1019. etarldouble (r, l)
  1020.      REAL_VALUE_TYPE r;
  1021.      long l[];
  1022. {
  1023.   unsigned EMUSHORT e[NE];
  1024.  
  1025.   GET_REAL (&r, e);
  1026.   etoe64 (e, e);
  1027.   endian (e, l, XFmode);
  1028. }
  1029.  
  1030. void 
  1031. etardouble (r, l)
  1032.      REAL_VALUE_TYPE r;
  1033.      long l[];
  1034. {
  1035.   unsigned EMUSHORT e[NE];
  1036.  
  1037.   GET_REAL (&r, e);
  1038.   etoe53 (e, e);
  1039.   endian (e, l, DFmode);
  1040. }
  1041.  
  1042. long
  1043. etarsingle (r)
  1044.      REAL_VALUE_TYPE r;
  1045. {
  1046.   unsigned EMUSHORT e[NE];
  1047.   long l;
  1048.  
  1049.   GET_REAL (&r, e);
  1050.   etoe24 (e, e);
  1051.   endian (e, &l, SFmode);
  1052.   return ((long) l);
  1053. }
  1054.  
  1055. void
  1056. ereal_to_decimal (x, s)
  1057.      REAL_VALUE_TYPE x;
  1058.      char *s;
  1059. {
  1060.   unsigned EMUSHORT e[NE];
  1061.  
  1062.   GET_REAL (&x, e);
  1063.   etoasc (e, s, 20);
  1064. }
  1065.  
  1066. int
  1067. ereal_cmp (x, y)
  1068.      REAL_VALUE_TYPE x, y;
  1069. {
  1070.   unsigned EMUSHORT ex[NE], ey[NE];
  1071.  
  1072.   GET_REAL (&x, ex);
  1073.   GET_REAL (&y, ey);
  1074.   return (ecmp (ex, ey));
  1075. }
  1076.  
  1077. int
  1078. ereal_isneg (x)
  1079.      REAL_VALUE_TYPE x;
  1080. {
  1081.   unsigned EMUSHORT ex[NE];
  1082.  
  1083.   GET_REAL (&x, ex);
  1084.   return (eisneg (ex));
  1085. }
  1086.  
  1087. /* End of REAL_ARITHMETIC interface */
  1088.  
  1089. /*
  1090.   Extended precision IEEE binary floating point arithmetic routines
  1091.  
  1092.   Numbers are stored in C language as arrays of 16-bit unsigned
  1093.   short integers.  The arguments of the routines are pointers to
  1094.   the arrays.
  1095.  
  1096.   External e type data structure, simulates Intel 8087 chip
  1097.   temporary real format but possibly with a larger significand:
  1098.  
  1099.     NE-1 significand words    (least significant word first,
  1100.                  most significant bit is normally set)
  1101.     exponent        (value = EXONE for 1.0,
  1102.                 top bit is the sign)
  1103.  
  1104.  
  1105.   Internal data structure of a number (a "word" is 16 bits):
  1106.  
  1107.   ei[0]    sign word    (0 for positive, 0xffff for negative)
  1108.   ei[1]    biased exponent    (value = EXONE for the number 1.0)
  1109.   ei[2]    high guard word    (always zero after normalization)
  1110.   ei[3]
  1111.   to ei[NI-2]    significand    (NI-4 significand words,
  1112.                   most significant word first,
  1113.                   most significant bit is set)
  1114.   ei[NI-1]    low guard word    (0x8000 bit is rounding place)
  1115.  
  1116.  
  1117.  
  1118.          Routines for external format numbers
  1119.  
  1120.      asctoe (string, e)    ASCII string to extended double e type
  1121.      asctoe64 (string, &d)    ASCII string to long double
  1122.      asctoe53 (string, &d)    ASCII string to double
  1123.      asctoe24 (string, &f)    ASCII string to single
  1124.      asctoeg (string, e, prec) ASCII string to specified precision
  1125.      e24toe (&f, e)        IEEE single precision to e type
  1126.      e53toe (&d, e)        IEEE double precision to e type
  1127.      e64toe (&d, e)        IEEE long double precision to e type
  1128.      e113toe (&d, e)        128-bit long double precision to e type
  1129.      eabs (e)            absolute value
  1130.      eadd (a, b, c)        c = b + a
  1131.      eclear (e)        e = 0
  1132.      ecmp (a, b)        Returns 1 if a > b, 0 if a == b,
  1133.                  -1 if a < b, -2 if either a or b is a NaN.
  1134.      ediv (a, b, c)        c = b / a
  1135.      efloor (a, b)        truncate to integer, toward -infinity
  1136.      efrexp (a, exp, s)    extract exponent and significand
  1137.      eifrac (e, &l, frac)    e to HOST_WIDE_INT and e type fraction
  1138.      euifrac (e, &l, frac)   e to unsigned HOST_WIDE_INT and e type fraction
  1139.      einfin (e)        set e to infinity, leaving its sign alone
  1140.      eldexp (a, n, b)    multiply by 2**n
  1141.      emov (a, b)        b = a
  1142.      emul (a, b, c)        c = b * a
  1143.      eneg (e)            e = -e
  1144.      eround (a, b)        b = nearest integer value to a
  1145.      esub (a, b, c)        c = b - a
  1146.      e24toasc (&f, str, n)    single to ASCII string, n digits after decimal
  1147.      e53toasc (&d, str, n)    double to ASCII string, n digits after decimal
  1148.      e64toasc (&d, str, n)    80-bit long double to ASCII string
  1149.      e113toasc (&d, str, n)    128-bit long double to ASCII string
  1150.      etoasc (e, str, n)    e to ASCII string, n digits after decimal
  1151.      etoe24 (e, &f)        convert e type to IEEE single precision
  1152.      etoe53 (e, &d)        convert e type to IEEE double precision
  1153.      etoe64 (e, &d)        convert e type to IEEE long double precision
  1154.      ltoe (&l, e)        HOST_WIDE_INT to e type
  1155.      ultoe (&l, e)        unsigned HOST_WIDE_INT to e type
  1156.     eisneg (e)              1 if sign bit of e != 0, else 0
  1157.     eisinf (e)              1 if e has maximum exponent (non-IEEE)
  1158.                  or is infinite (IEEE)
  1159.         eisnan (e)              1 if e is a NaN
  1160.  
  1161.  
  1162.          Routines for internal format numbers
  1163.  
  1164.      eaddm (ai, bi)        add significands, bi = bi + ai
  1165.      ecleaz (ei)        ei = 0
  1166.      ecleazs (ei)        set ei = 0 but leave its sign alone
  1167.      ecmpm (ai, bi)        compare significands, return 1, 0, or -1
  1168.      edivm (ai, bi)        divide  significands, bi = bi / ai
  1169.      emdnorm (ai,l,s,exp)    normalize and round off
  1170.      emovi (a, ai)        convert external a to internal ai
  1171.      emovo (ai, a)        convert internal ai to external a
  1172.      emovz (ai, bi)        bi = ai, low guard word of bi = 0
  1173.      emulm (ai, bi)        multiply significands, bi = bi * ai
  1174.      enormlz (ei)        left-justify the significand
  1175.      eshdn1 (ai)        shift significand and guards down 1 bit
  1176.      eshdn8 (ai)        shift down 8 bits
  1177.      eshdn6 (ai)        shift down 16 bits
  1178.      eshift (ai, n)        shift ai n bits up (or down if n < 0)
  1179.      eshup1 (ai)        shift significand and guards up 1 bit
  1180.      eshup8 (ai)        shift up 8 bits
  1181.      eshup6 (ai)        shift up 16 bits
  1182.      esubm (ai, bi)        subtract significands, bi = bi - ai
  1183.         eiisinf (ai)            1 if infinite
  1184.         eiisnan (ai)            1 if a NaN
  1185.      eiisneg (ai)        1 if sign bit of ai != 0, else 0
  1186.         einan (ai)              set ai = NaN
  1187.         eiinfin (ai)            set ai = infinity
  1188.  
  1189.   The result is always normalized and rounded to NI-4 word precision
  1190.   after each arithmetic operation.
  1191.  
  1192.   Exception flags are NOT fully supported.
  1193.  
  1194.   Signaling NaN's are NOT supported; they are treated the same
  1195.   as quiet NaN's.
  1196.  
  1197.   Define INFINITY for support of infinity; otherwise a
  1198.   saturation arithmetic is implemented.
  1199.  
  1200.   Define NANS for support of Not-a-Number items; otherwise the
  1201.   arithmetic will never produce a NaN output, and might be confused
  1202.   by a NaN input.
  1203.   If NaN's are supported, the output of `ecmp (a,b)' is -2 if
  1204.   either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
  1205.   may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
  1206.   if in doubt.
  1207.  
  1208.   Denormals are always supported here where appropriate (e.g., not
  1209.   for conversion to DEC numbers).  */
  1210.  
  1211. /* Definitions for error codes that are passed to the common error handling
  1212.    routine mtherr.
  1213.  
  1214.    For Digital Equipment PDP-11 and VAX computers, certain
  1215.   IBM systems, and others that use numbers with a 56-bit
  1216.   significand, the symbol DEC should be defined.  In this
  1217.   mode, most floating point constants are given as arrays
  1218.   of octal integers to eliminate decimal to binary conversion
  1219.   errors that might be introduced by the compiler.
  1220.  
  1221.   For computers, such as IBM PC, that follow the IEEE
  1222.   Standard for Binary Floating Point Arithmetic (ANSI/IEEE
  1223.   Std 754-1985), the symbol IBMPC or MIEEE should be defined.
  1224.   These numbers have 53-bit significands.  In this mode, constants
  1225.   are provided as arrays of hexadecimal 16 bit integers.
  1226.  
  1227.   To accommodate other types of computer arithmetic, all
  1228.   constants are also provided in a normal decimal radix
  1229.   which one can hope are correctly converted to a suitable
  1230.   format by the available C language compiler.  To invoke
  1231.   this mode, the symbol UNK is defined.
  1232.  
  1233.   An important difference among these modes is a predefined
  1234.   set of machine arithmetic constants for each.  The numbers
  1235.   MACHEP (the machine roundoff error), MAXNUM (largest number
  1236.   represented), and several other parameters are preset by
  1237.   the configuration symbol.  Check the file const.c to
  1238.   ensure that these values are correct for your computer.
  1239.  
  1240.   For ANSI C compatibility, define ANSIC equal to 1.  Currently
  1241.   this affects only the atan2 function and others that use it. */
  1242.  
  1243. /* Constant definitions for math error conditions.  */
  1244.  
  1245. #define DOMAIN        1    /* argument domain error */
  1246. #define SING        2    /* argument singularity */
  1247. #define OVERFLOW    3    /* overflow range error */
  1248. #define UNDERFLOW    4    /* underflow range error */
  1249. #define TLOSS        5    /* total loss of precision */
  1250. #define PLOSS        6    /* partial loss of precision */
  1251. #define INVALID        7    /* NaN-producing operation */
  1252.  
  1253. /*  e type constants used by high precision check routines */
  1254.  
  1255. #if LONG_DOUBLE_TYPE_SIZE == 128
  1256. /* 0.0 */
  1257. unsigned EMUSHORT ezero[NE] =
  1258.  {0x0000, 0x0000, 0x0000, 0x0000,
  1259.   0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
  1260. extern unsigned EMUSHORT ezero[];
  1261.  
  1262. /* 5.0E-1 */
  1263. unsigned EMUSHORT ehalf[NE] =
  1264.  {0x0000, 0x0000, 0x0000, 0x0000,
  1265.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
  1266. extern unsigned EMUSHORT ehalf[];
  1267.  
  1268. /* 1.0E0 */
  1269. unsigned EMUSHORT eone[NE] =
  1270.  {0x0000, 0x0000, 0x0000, 0x0000,
  1271.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
  1272. extern unsigned EMUSHORT eone[];
  1273.  
  1274. /* 2.0E0 */
  1275. unsigned EMUSHORT etwo[NE] =
  1276.  {0x0000, 0x0000, 0x0000, 0x0000,
  1277.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
  1278. extern unsigned EMUSHORT etwo[];
  1279.  
  1280. /* 3.2E1 */
  1281. unsigned EMUSHORT e32[NE] =
  1282.  {0x0000, 0x0000, 0x0000, 0x0000,
  1283.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
  1284. extern unsigned EMUSHORT e32[];
  1285.  
  1286. /* 6.93147180559945309417232121458176568075500134360255E-1 */
  1287. unsigned EMUSHORT elog2[NE] =
  1288.  {0x40f3, 0xf6af, 0x03f2, 0xb398,
  1289.   0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
  1290. extern unsigned EMUSHORT elog2[];
  1291.  
  1292. /* 1.41421356237309504880168872420969807856967187537695E0 */
  1293. unsigned EMUSHORT esqrt2[NE] =
  1294.  {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
  1295.   0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
  1296. extern unsigned EMUSHORT esqrt2[];
  1297.  
  1298. /* 3.14159265358979323846264338327950288419716939937511E0 */
  1299. unsigned EMUSHORT epi[NE] =
  1300.  {0x2902, 0x1cd1, 0x80dc, 0x628b,
  1301.   0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
  1302. extern unsigned EMUSHORT epi[];
  1303.  
  1304. #else
  1305. /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
  1306. unsigned EMUSHORT ezero[NE] =
  1307.  {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
  1308. unsigned EMUSHORT ehalf[NE] =
  1309.  {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
  1310. unsigned EMUSHORT eone[NE] =
  1311.  {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
  1312. unsigned EMUSHORT etwo[NE] =
  1313.  {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
  1314. unsigned EMUSHORT e32[NE] =
  1315.  {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
  1316. unsigned EMUSHORT elog2[NE] =
  1317.  {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
  1318. unsigned EMUSHORT esqrt2[NE] =
  1319.  {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
  1320. unsigned EMUSHORT epi[NE] =
  1321.  {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
  1322. #endif
  1323.  
  1324.  
  1325.  
  1326. /* Control register for rounding precision.
  1327.    This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits.  */
  1328.  
  1329. int rndprc = NBITS;
  1330. extern int rndprc;
  1331.  
  1332. /*  Clear out entire external format number.  */
  1333.  
  1334. static void 
  1335. eclear (x)
  1336.      register unsigned EMUSHORT *x;
  1337. {
  1338.   register int i;
  1339.  
  1340.   for (i = 0; i < NE; i++)
  1341.     *x++ = 0;
  1342. }
  1343.  
  1344.  
  1345.  
  1346. /* Move external format number from a to b.  */
  1347.  
  1348. static void 
  1349. emov (a, b)
  1350.      register unsigned EMUSHORT *a, *b;
  1351. {
  1352.   register int i;
  1353.  
  1354.   for (i = 0; i < NE; i++)
  1355.     *b++ = *a++;
  1356. }
  1357.  
  1358.  
  1359. /* Absolute value of external format number.  */
  1360.  
  1361. static void 
  1362. eabs (x)
  1363.      unsigned EMUSHORT x[];
  1364. {
  1365.   /* sign is top bit of last word of external format */
  1366.   x[NE - 1] &= 0x7fff;        
  1367. }
  1368.  
  1369. /* Negate external format number.  */
  1370.  
  1371. static void 
  1372. eneg (x)
  1373.      unsigned EMUSHORT x[];
  1374. {
  1375.  
  1376.   x[NE - 1] ^= 0x8000;        /* Toggle the sign bit */
  1377. }
  1378.  
  1379.  
  1380.  
  1381. /* Return 1 if sign bit of external format number is nonzero, else zero.  */
  1382.  
  1383. static int 
  1384. eisneg (x)
  1385.      unsigned EMUSHORT x[];
  1386. {
  1387.  
  1388.   if (x[NE - 1] & 0x8000)
  1389.     return (1);
  1390.   else
  1391.     return (0);
  1392. }
  1393.  
  1394.  
  1395. /* Return 1 if external format number is infinity, else return zero.  */
  1396.  
  1397. static int 
  1398. eisinf (x)
  1399.      unsigned EMUSHORT x[];
  1400. {
  1401.  
  1402. #ifdef NANS
  1403.   if (eisnan (x))
  1404.     return (0);
  1405. #endif
  1406.   if ((x[NE - 1] & 0x7fff) == 0x7fff)
  1407.     return (1);
  1408.   else
  1409.     return (0);
  1410. }
  1411.  
  1412.  
  1413. /* Check if e-type number is not a number.  The bit pattern is one that we
  1414.    defined, so we know for sure how to detect it.  */
  1415.  
  1416. static int 
  1417. eisnan (x)
  1418.      unsigned EMUSHORT x[];
  1419. {
  1420. #ifdef NANS
  1421.   int i;
  1422.  
  1423.   /* NaN has maximum exponent */
  1424.   if ((x[NE - 1] & 0x7fff) != 0x7fff)
  1425.     return (0);
  1426.   /* ... and non-zero significand field. */
  1427.   for (i = 0; i < NE - 1; i++)
  1428.     {
  1429.       if (*x++ != 0)
  1430.         return (1);
  1431.     }
  1432. #endif
  1433.  
  1434.   return (0);
  1435. }
  1436.  
  1437. /*  Fill external format number with infinity pattern (IEEE)
  1438.     or largest possible number (non-IEEE). */
  1439.  
  1440. static void 
  1441. einfin (x)
  1442.      register unsigned EMUSHORT *x;
  1443. {
  1444.   register int i;
  1445.  
  1446. #ifdef INFINITY
  1447.   for (i = 0; i < NE - 1; i++)
  1448.     *x++ = 0;
  1449.   *x |= 32767;
  1450. #else
  1451.   for (i = 0; i < NE - 1; i++)
  1452.     *x++ = 0xffff;
  1453.   *x |= 32766;
  1454.   if (rndprc < NBITS)
  1455.     {
  1456.       if (rndprc == 113)
  1457.     {
  1458.       *(x - 9) = 0;
  1459.       *(x - 8) = 0;
  1460.     }
  1461.       if (rndprc == 64)
  1462.     {
  1463.       *(x - 5) = 0;
  1464.     }
  1465.       if (rndprc == 53)
  1466.     {
  1467.       *(x - 4) = 0xf800;
  1468.     }
  1469.       else
  1470.     {
  1471.       *(x - 4) = 0;
  1472.       *(x - 3) = 0;
  1473.       *(x - 2) = 0xff00;
  1474.     }
  1475.     }
  1476. #endif
  1477. }
  1478.  
  1479.  
  1480. /* Output an e-type NaN.
  1481.    This generates Intel's quiet NaN pattern for extended real.
  1482.    The exponent is 7fff, the leading mantissa word is c000.  */
  1483.  
  1484. static void 
  1485. enan (x, sign)
  1486.      register unsigned EMUSHORT *x;
  1487.      int sign;
  1488. {
  1489.   register int i;
  1490.  
  1491.   for (i = 0; i < NE - 2; i++)
  1492.     *x++ = 0;
  1493.   *x++ = 0xc000;
  1494.   *x = (sign << 15) | 0x7fff;
  1495. }
  1496.  
  1497.  
  1498. /* Move in external format number, converting it to internal format.  */
  1499.  
  1500. static void 
  1501. emovi (a, b)
  1502.      unsigned EMUSHORT *a, *b;
  1503. {
  1504.   register unsigned EMUSHORT *p, *q;
  1505.   int i;
  1506.  
  1507.   q = b;
  1508.   p = a + (NE - 1);        /* point to last word of external number */
  1509.   /* get the sign bit */
  1510.   if (*p & 0x8000)
  1511.     *q++ = 0xffff;
  1512.   else
  1513.     *q++ = 0;
  1514.   /* get the exponent */
  1515.   *q = *p--;
  1516.   *q++ &= 0x7fff;        /* delete the sign bit */
  1517. #ifdef INFINITY
  1518.   if ((*(q - 1) & 0x7fff) == 0x7fff)
  1519.     {
  1520. #ifdef NANS
  1521.       if (eisnan (a))
  1522.     {
  1523.       *q++ = 0;
  1524.       for (i = 3; i < NI; i++)
  1525.         *q++ = *p--;
  1526.       return;
  1527.     }
  1528. #endif
  1529.  
  1530.       for (i = 2; i < NI; i++)
  1531.     *q++ = 0;
  1532.       return;
  1533.     }
  1534. #endif
  1535.  
  1536.   /* clear high guard word */
  1537.   *q++ = 0;
  1538.   /* move in the significand */
  1539.   for (i = 0; i < NE - 1; i++)
  1540.     *q++ = *p--;
  1541.   /* clear low guard word */
  1542.   *q = 0;
  1543. }
  1544.  
  1545.  
  1546. /* Move internal format number out, converting it to external format.  */
  1547.  
  1548. static void 
  1549. emovo (a, b)
  1550.      unsigned EMUSHORT *a, *b;
  1551. {
  1552.   register unsigned EMUSHORT *p, *q;
  1553.   unsigned EMUSHORT i;
  1554.   int j;
  1555.  
  1556.   p = a;
  1557.   q = b + (NE - 1);        /* point to output exponent */
  1558.   /* combine sign and exponent */
  1559.   i = *p++;
  1560.   if (i)
  1561.     *q-- = *p++ | 0x8000;
  1562.   else
  1563.     *q-- = *p++;
  1564. #ifdef INFINITY
  1565.   if (*(p - 1) == 0x7fff)
  1566.     {
  1567. #ifdef NANS
  1568.       if (eiisnan (a))
  1569.     {
  1570.       enan (b, eiisneg (a));
  1571.       return;
  1572.     }
  1573. #endif
  1574.       einfin (b);
  1575.     return;
  1576.     }
  1577. #endif
  1578.   /* skip over guard word */
  1579.   ++p;
  1580.   /* move the significand */
  1581.   for (j = 0; j < NE - 1; j++)
  1582.     *q-- = *p++;
  1583. }
  1584.  
  1585. /* Clear out internal format number.  */
  1586.  
  1587. static void 
  1588. ecleaz (xi)
  1589.      register unsigned EMUSHORT *xi;
  1590. {
  1591.   register int i;
  1592.  
  1593.   for (i = 0; i < NI; i++)
  1594.     *xi++ = 0;
  1595. }
  1596.  
  1597.  
  1598. /* Same, but don't touch the sign. */
  1599.  
  1600. static void 
  1601. ecleazs (xi)
  1602.      register unsigned EMUSHORT *xi;
  1603. {
  1604.   register int i;
  1605.  
  1606.   ++xi;
  1607.   for (i = 0; i < NI - 1; i++)
  1608.     *xi++ = 0;
  1609. }
  1610.  
  1611.  
  1612.  
  1613. /* Move internal format number from a to b.  */
  1614.  
  1615. static void 
  1616. emovz (a, b)
  1617.      register unsigned EMUSHORT *a, *b;
  1618. {
  1619.   register int i;
  1620.  
  1621.   for (i = 0; i < NI - 1; i++)
  1622.     *b++ = *a++;
  1623.   /* clear low guard word */
  1624.   *b = 0;
  1625. }
  1626.  
  1627. /* Generate internal format NaN.
  1628.    The explicit pattern for this is maximum exponent and
  1629.    top two significant bits set.  */
  1630.  
  1631. static void
  1632. einan (x)
  1633.      unsigned EMUSHORT x[];
  1634. {
  1635.  
  1636.   ecleaz (x);
  1637.   x[E] = 0x7fff;
  1638.   x[M + 1] = 0xc000;
  1639. }
  1640.  
  1641. /* Return nonzero if internal format number is a NaN. */
  1642.  
  1643. static int 
  1644. eiisnan (x)
  1645.      unsigned EMUSHORT x[];
  1646. {
  1647.   int i;
  1648.  
  1649.   if ((x[E] & 0x7fff) == 0x7fff)
  1650.     {
  1651.       for (i = M + 1; i < NI; i++)
  1652.     {
  1653.       if (x[i] != 0)
  1654.         return (1);
  1655.     }
  1656.     }
  1657.   return (0);
  1658. }
  1659.  
  1660. /* Return nonzero if sign of internal format number is nonzero.  */
  1661.  
  1662. static int 
  1663. eiisneg (x)
  1664.      unsigned EMUSHORT x[];
  1665. {
  1666.  
  1667.   return x[0] != 0;
  1668. }
  1669.  
  1670. /* Fill internal format number with infinity pattern.
  1671.    This has maximum exponent and significand all zeros.  */
  1672.  
  1673. static void
  1674. eiinfin (x)
  1675.      unsigned EMUSHORT x[];
  1676. {
  1677.  
  1678.   ecleaz (x);
  1679.   x[E] = 0x7fff;
  1680. }
  1681.  
  1682. /* Return nonzero if internal format number is infinite. */
  1683.  
  1684. static int 
  1685. eiisinf (x)
  1686.      unsigned EMUSHORT x[];
  1687. {
  1688.  
  1689. #ifdef NANS
  1690.   if (eiisnan (x))
  1691.     return (0);
  1692. #endif
  1693.   if ((x[E] & 0x7fff) == 0x7fff)
  1694.     return (1);
  1695.   return (0);
  1696. }
  1697.  
  1698.  
  1699. /* Compare significands of numbers in internal format.
  1700.    Guard words are included in the comparison.
  1701.  
  1702.    Returns    +1 if a > b
  1703.          0 if a == b
  1704.         -1 if a < b   */
  1705.  
  1706. static int
  1707. ecmpm (a, b)
  1708.      register unsigned EMUSHORT *a, *b;
  1709. {
  1710.   int i;
  1711.  
  1712.   a += M;            /* skip up to significand area */
  1713.   b += M;
  1714.   for (i = M; i < NI; i++)
  1715.     {
  1716.       if (*a++ != *b++)
  1717.     goto difrnt;
  1718.     }
  1719.   return (0);
  1720.  
  1721.  difrnt:
  1722.   if (*(--a) > *(--b))
  1723.     return (1);
  1724.   else
  1725.     return (-1);
  1726. }
  1727.  
  1728.  
  1729. /* Shift significand down by 1 bit.  */
  1730.  
  1731. static void 
  1732. eshdn1 (x)
  1733.      register unsigned EMUSHORT *x;
  1734. {
  1735.   register unsigned EMUSHORT bits;
  1736.   int i;
  1737.  
  1738.   x += M;            /* point to significand area */
  1739.  
  1740.   bits = 0;
  1741.   for (i = M; i < NI; i++)
  1742.     {
  1743.       if (*x & 1)
  1744.     bits |= 1;
  1745.       *x >>= 1;
  1746.       if (bits & 2)
  1747.     *x |= 0x8000;
  1748.       bits <<= 1;
  1749.       ++x;
  1750.     }
  1751. }
  1752.  
  1753.  
  1754.  
  1755. /* Shift significand up by 1 bit.  */
  1756.  
  1757. static void 
  1758. eshup1 (x)
  1759.      register unsigned EMUSHORT *x;
  1760. {
  1761.   register unsigned EMUSHORT bits;
  1762.   int i;
  1763.  
  1764.   x += NI - 1;
  1765.   bits = 0;
  1766.  
  1767.   for (i = M; i < NI; i++)
  1768.     {
  1769.       if (*x & 0x8000)
  1770.     bits |= 1;
  1771.       *x <<= 1;
  1772.       if (bits & 2)
  1773.     *x |= 1;
  1774.       bits <<= 1;
  1775.       --x;
  1776.     }
  1777. }
  1778.  
  1779.  
  1780. /* Shift significand down by 8 bits.  */
  1781.  
  1782. static void 
  1783. eshdn8 (x)
  1784.      register unsigned EMUSHORT *x;
  1785. {
  1786.   register unsigned EMUSHORT newbyt, oldbyt;
  1787.   int i;
  1788.  
  1789.   x += M;
  1790.   oldbyt = 0;
  1791.   for (i = M; i < NI; i++)
  1792.     {
  1793.       newbyt = *x << 8;
  1794.       *x >>= 8;
  1795.       *x |= oldbyt;
  1796.       oldbyt = newbyt;
  1797.       ++x;
  1798.     }
  1799. }
  1800.  
  1801. /* Shift significand up by 8 bits.  */
  1802.  
  1803. static void 
  1804. eshup8 (x)
  1805.      register unsigned EMUSHORT *x;
  1806. {
  1807.   int i;
  1808.   register unsigned EMUSHORT newbyt, oldbyt;
  1809.  
  1810.   x += NI - 1;
  1811.   oldbyt = 0;
  1812.  
  1813.   for (i = M; i < NI; i++)
  1814.     {
  1815.       newbyt = *x >> 8;
  1816.       *x <<= 8;
  1817.       *x |= oldbyt;
  1818.       oldbyt = newbyt;
  1819.       --x;
  1820.     }
  1821. }
  1822.  
  1823. /* Shift significand up by 16 bits.  */
  1824.  
  1825. static void 
  1826. eshup6 (x)
  1827.      register unsigned EMUSHORT *x;
  1828. {
  1829.   int i;
  1830.   register unsigned EMUSHORT *p;
  1831.  
  1832.   p = x + M;
  1833.   x += M + 1;
  1834.  
  1835.   for (i = M; i < NI - 1; i++)
  1836.     *p++ = *x++;
  1837.  
  1838.   *p = 0;
  1839. }
  1840.  
  1841. /* Shift significand down by 16 bits.  */
  1842.  
  1843. static void 
  1844. eshdn6 (x)
  1845.      register unsigned EMUSHORT *x;
  1846. {
  1847.   int i;
  1848.   register unsigned EMUSHORT *p;
  1849.  
  1850.   x += NI - 1;
  1851.   p = x + 1;
  1852.  
  1853.   for (i = M; i < NI - 1; i++)
  1854.     *(--p) = *(--x);
  1855.  
  1856.   *(--p) = 0;
  1857. }
  1858.  
  1859. /* Add significands.  x + y replaces y.  */
  1860.  
  1861. static void 
  1862. eaddm (x, y)
  1863.      unsigned EMUSHORT *x, *y;
  1864. {
  1865.   register unsigned EMULONG a;
  1866.   int i;
  1867.   unsigned int carry;
  1868.  
  1869.   x += NI - 1;
  1870.   y += NI - 1;
  1871.   carry = 0;
  1872.   for (i = M; i < NI; i++)
  1873.     {
  1874.       a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
  1875.       if (a & 0x10000)
  1876.     carry = 1;
  1877.       else
  1878.     carry = 0;
  1879.       *y = (unsigned EMUSHORT) a;
  1880.       --x;
  1881.       --y;
  1882.     }
  1883. }
  1884.  
  1885. /* Subtract significands.  y - x replaces y.  */
  1886.  
  1887. static void 
  1888. esubm (x, y)
  1889.      unsigned EMUSHORT *x, *y;
  1890. {
  1891.   unsigned EMULONG a;
  1892.   int i;
  1893.   unsigned int carry;
  1894.  
  1895.   x += NI - 1;
  1896.   y += NI - 1;
  1897.   carry = 0;
  1898.   for (i = M; i < NI; i++)
  1899.     {
  1900.       a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
  1901.       if (a & 0x10000)
  1902.     carry = 1;
  1903.       else
  1904.     carry = 0;
  1905.       *y = (unsigned EMUSHORT) a;
  1906.       --x;
  1907.       --y;
  1908.     }
  1909. }
  1910.  
  1911.  
  1912. static unsigned EMUSHORT equot[NI];
  1913.  
  1914.  
  1915. #if 0
  1916. /* Radix 2 shift-and-add versions of multiply and divide  */
  1917.  
  1918.  
  1919. /* Divide significands */
  1920.  
  1921. int 
  1922. edivm (den, num)
  1923.      unsigned EMUSHORT den[], num[];
  1924. {
  1925.   int i;
  1926.   register unsigned EMUSHORT *p, *q;
  1927.   unsigned EMUSHORT j;
  1928.  
  1929.   p = &equot[0];
  1930.   *p++ = num[0];
  1931.   *p++ = num[1];
  1932.  
  1933.   for (i = M; i < NI; i++)
  1934.     {
  1935.       *p++ = 0;
  1936.     }
  1937.  
  1938.   /* Use faster compare and subtraction if denominator has only 15 bits of
  1939.      significance.  */
  1940.  
  1941.   p = &den[M + 2];
  1942.   if (*p++ == 0)
  1943.     {
  1944.       for (i = M + 3; i < NI; i++)
  1945.     {
  1946.       if (*p++ != 0)
  1947.         goto fulldiv;
  1948.     }
  1949.       if ((den[M + 1] & 1) != 0)
  1950.     goto fulldiv;
  1951.       eshdn1 (num);
  1952.       eshdn1 (den);
  1953.  
  1954.       p = &den[M + 1];
  1955.       q = &num[M + 1];
  1956.  
  1957.       for (i = 0; i < NBITS + 2; i++)
  1958.     {
  1959.       if (*p <= *q)
  1960.         {
  1961.           *q -= *p;
  1962.           j = 1;
  1963.         }
  1964.       else
  1965.         {
  1966.           j = 0;
  1967.         }
  1968.       eshup1 (equot);
  1969.       equot[NI - 2] |= j;
  1970.       eshup1 (num);
  1971.     }
  1972.       goto divdon;
  1973.     }
  1974.  
  1975.   /* The number of quotient bits to calculate is NBITS + 1 scaling guard
  1976.      bit + 1 roundoff bit.  */
  1977.  
  1978.  fulldiv:
  1979.  
  1980.   p = &equot[NI - 2];
  1981.   for (i = 0; i < NBITS + 2; i++)
  1982.     {
  1983.       if (ecmpm (den, num) <= 0)
  1984.     {
  1985.       esubm (den, num);
  1986.       j = 1;        /* quotient bit = 1 */
  1987.     }
  1988.       else
  1989.     j = 0;
  1990.       eshup1 (equot);
  1991.       *p |= j;
  1992.       eshup1 (num);
  1993.     }
  1994.  
  1995.  divdon:
  1996.  
  1997.   eshdn1 (equot);
  1998.   eshdn1 (equot);
  1999.  
  2000.   /* test for nonzero remainder after roundoff bit */
  2001.   p = &num[M];
  2002.   j = 0;
  2003.   for (i = M; i < NI; i++)
  2004.     {
  2005.       j |= *p++;
  2006.     }
  2007.   if (j)
  2008.     j = 1;
  2009.  
  2010.  
  2011.   for (i = 0; i < NI; i++)
  2012.     num[i] = equot[i];
  2013.   return ((int) j);
  2014. }
  2015.  
  2016.  
  2017. /* Multiply significands */
  2018. int 
  2019. emulm (a, b)
  2020.      unsigned EMUSHORT a[], b[];
  2021. {
  2022.   unsigned EMUSHORT *p, *q;
  2023.   int i, j, k;
  2024.  
  2025.   equot[0] = b[0];
  2026.   equot[1] = b[1];
  2027.   for (i = M; i < NI; i++)
  2028.     equot[i] = 0;
  2029.  
  2030.   p = &a[NI - 2];
  2031.   k = NBITS;
  2032.   while (*p == 0)        /* significand is not supposed to be zero */
  2033.     {
  2034.       eshdn6 (a);
  2035.       k -= 16;
  2036.     }
  2037.   if ((*p & 0xff) == 0)
  2038.     {
  2039.       eshdn8 (a);
  2040.       k -= 8;
  2041.     }
  2042.  
  2043.   q = &equot[NI - 1];
  2044.   j = 0;
  2045.   for (i = 0; i < k; i++)
  2046.     {
  2047.       if (*p & 1)
  2048.     eaddm (b, equot);
  2049.       /* remember if there were any nonzero bits shifted out */
  2050.       if (*q & 1)
  2051.     j |= 1;
  2052.       eshdn1 (a);
  2053.       eshdn1 (equot);
  2054.     }
  2055.  
  2056.   for (i = 0; i < NI; i++)
  2057.     b[i] = equot[i];
  2058.  
  2059.   /* return flag for lost nonzero bits */
  2060.   return (j);
  2061. }
  2062.  
  2063. #else
  2064.  
  2065. /* Radix 65536 versions of multiply and divide  */
  2066.  
  2067.  
  2068. /* Multiply significand of e-type number b
  2069.    by 16-bit quantity a, e-type result to c. */
  2070.  
  2071. static void
  2072. m16m (a, b, c)
  2073.      unsigned int a;
  2074.      unsigned short b[], c[];
  2075. {
  2076.   register unsigned short *pp;
  2077.   register unsigned long carry;
  2078.   unsigned short *ps;
  2079.   unsigned short p[NI];
  2080.   unsigned long aa, m;
  2081.   int i;
  2082.  
  2083.   aa = a;
  2084.   pp = &p[NI-2];
  2085.   *pp++ = 0;
  2086.   *pp = 0;
  2087.   ps = &b[NI-1];
  2088.  
  2089.   for (i=M+1; i<NI; i++)
  2090.     {
  2091.       if (*ps == 0)
  2092.     {
  2093.       --ps;
  2094.       --pp;
  2095.       *(pp-1) = 0;
  2096.     }
  2097.       else
  2098.     {
  2099.       m = (unsigned long) aa * *ps--;
  2100.       carry = (m & 0xffff) + *pp;
  2101.       *pp-- = (unsigned short)carry;
  2102.       carry = (carry >> 16) + (m >> 16) + *pp;
  2103.       *pp = (unsigned short)carry;
  2104.       *(pp-1) = carry >> 16;
  2105.     }
  2106.     }
  2107.   for (i=M; i<NI; i++)
  2108.     c[i] = p[i];
  2109. }
  2110.  
  2111.  
  2112. /* Divide significands. Neither the numerator nor the denominator
  2113.    is permitted to have its high guard word nonzero.  */
  2114.  
  2115. static int
  2116. edivm (den, num)
  2117.      unsigned short den[], num[];
  2118. {
  2119.   int i;
  2120.   register unsigned short *p;
  2121.   unsigned long tnum;
  2122.   unsigned short j, tdenm, tquot;
  2123.   unsigned short tprod[NI+1];
  2124.  
  2125.   p = &equot[0];
  2126.   *p++ = num[0];
  2127.   *p++ = num[1];
  2128.  
  2129.   for (i=M; i<NI; i++)
  2130.     {
  2131.       *p++ = 0;
  2132.     }
  2133.   eshdn1 (num);
  2134.   tdenm = den[M+1];
  2135.   for (i=M; i<NI; i++)
  2136.     {
  2137.       /* Find trial quotient digit (the radix is 65536). */
  2138.       tnum = (((unsigned long) num[M]) << 16) + num[M+1];
  2139.  
  2140.       /* Do not execute the divide instruction if it will overflow. */
  2141.       if ((tdenm * 0xffffL) < tnum)
  2142.     tquot = 0xffff;
  2143.       else
  2144.     tquot = tnum / tdenm;
  2145.       /* Multiply denominator by trial quotient digit. */
  2146.       m16m ((unsigned int)tquot, den, tprod);
  2147.       /* The quotient digit may have been overestimated. */
  2148.       if (ecmpm (tprod, num) > 0)
  2149.     {
  2150.       tquot -= 1;
  2151.       esubm (den, tprod);
  2152.       if (ecmpm (tprod, num) > 0)
  2153.         {
  2154.           tquot -= 1;
  2155.           esubm (den, tprod);
  2156.         }
  2157.     }
  2158.       esubm (tprod, num);
  2159.       equot[i] = tquot;
  2160.       eshup6(num);
  2161.     }
  2162.   /* test for nonzero remainder after roundoff bit */
  2163.   p = &num[M];
  2164.   j = 0;
  2165.   for (i=M; i<NI; i++)
  2166.     {
  2167.       j |= *p++;
  2168.     }
  2169.   if (j)
  2170.     j = 1;
  2171.  
  2172.   for (i=0; i<NI; i++)
  2173.     num[i] = equot[i];
  2174.  
  2175.   return ((int)j);
  2176. }
  2177.  
  2178.  
  2179.  
  2180. /* Multiply significands */
  2181. static int
  2182. emulm (a, b)
  2183.      unsigned short a[], b[];
  2184. {
  2185.   unsigned short *p, *q;
  2186.   unsigned short pprod[NI];
  2187.   unsigned short j;
  2188.   int i;
  2189.  
  2190.   equot[0] = b[0];
  2191.   equot[1] = b[1];
  2192.   for (i=M; i<NI; i++)
  2193.     equot[i] = 0;
  2194.  
  2195.   j = 0;
  2196.   p = &a[NI-1];
  2197.   q = &equot[NI-1];
  2198.   for (i=M+1; i<NI; i++)
  2199.     {
  2200.       if (*p == 0)
  2201.     {
  2202.       --p;
  2203.     }
  2204.       else
  2205.     {
  2206.       m16m ((unsigned int) *p--, b, pprod);
  2207.       eaddm(pprod, equot);
  2208.     }
  2209.       j |= *q;
  2210.       eshdn6(equot);
  2211.     }
  2212.  
  2213.   for (i=0; i<NI; i++)
  2214.     b[i] = equot[i];
  2215.  
  2216.   /* return flag for lost nonzero bits */
  2217.   return ((int)j);
  2218. }
  2219. #endif
  2220.  
  2221.  
  2222. /* Normalize and round off.
  2223.  
  2224.   The internal format number to be rounded is "s".
  2225.   Input "lost" indicates whether or not the number is exact.
  2226.   This is the so-called sticky bit.
  2227.  
  2228.   Input "subflg" indicates whether the number was obtained
  2229.   by a subtraction operation.  In that case if lost is nonzero
  2230.   then the number is slightly smaller than indicated.
  2231.  
  2232.   Input "exp" is the biased exponent, which may be negative.
  2233.   the exponent field of "s" is ignored but is replaced by
  2234.   "exp" as adjusted by normalization and rounding.
  2235.  
  2236.   Input "rcntrl" is the rounding control.
  2237.  
  2238.   For future reference:  In order for emdnorm to round off denormal
  2239.    significands at the right point, the input exponent must be
  2240.    adjusted to be the actual value it would have after conversion to
  2241.    the final floating point type.  This adjustment has been
  2242.    implemented for all type conversions (etoe53, etc.) and decimal
  2243.    conversions, but not for the arithmetic functions (eadd, etc.). 
  2244.    Data types having standard 15-bit exponents are not affected by
  2245.    this, but SFmode and DFmode are affected. For example, ediv with
  2246.    rndprc = 24 will not round correctly to 24-bit precision if the
  2247.    result is denormal.   */
  2248.  
  2249. static int rlast = -1;
  2250. static int rw = 0;
  2251. static unsigned EMUSHORT rmsk = 0;
  2252. static unsigned EMUSHORT rmbit = 0;
  2253. static unsigned EMUSHORT rebit = 0;
  2254. static int re = 0;
  2255. static unsigned EMUSHORT rbit[NI];
  2256.  
  2257. static void 
  2258. emdnorm (s, lost, subflg, exp, rcntrl)
  2259.      unsigned EMUSHORT s[];
  2260.      int lost;
  2261.      int subflg;
  2262.      EMULONG exp;
  2263.      int rcntrl;
  2264. {
  2265.   int i, j;
  2266.   unsigned EMUSHORT r;
  2267.  
  2268.   /* Normalize */
  2269.   j = enormlz (s);
  2270.  
  2271.   /* a blank significand could mean either zero or infinity. */
  2272. #ifndef INFINITY
  2273.   if (j > NBITS)
  2274.     {
  2275.       ecleazs (s);
  2276.       return;
  2277.     }
  2278. #endif
  2279.   exp -= j;
  2280. #ifndef INFINITY
  2281.   if (exp >= 32767L)
  2282.     goto overf;
  2283. #else
  2284.   if ((j > NBITS) && (exp < 32767))
  2285.     {
  2286.       ecleazs (s);
  2287.       return;
  2288.     }
  2289. #endif
  2290.   if (exp < 0L)
  2291.     {
  2292.       if (exp > (EMULONG) (-NBITS - 1))
  2293.     {
  2294.       j = (int) exp;
  2295.       i = eshift (s, j);
  2296.       if (i)
  2297.         lost = 1;
  2298.     }
  2299.       else
  2300.     {
  2301.       ecleazs (s);
  2302.       return;
  2303.     }
  2304.     }
  2305.   /* Round off, unless told not to by rcntrl. */
  2306.   if (rcntrl == 0)
  2307.     goto mdfin;
  2308.   /* Set up rounding parameters if the control register changed. */
  2309.   if (rndprc != rlast)
  2310.     {
  2311.       ecleaz (rbit);
  2312.       switch (rndprc)
  2313.     {
  2314.     default:
  2315.     case NBITS:
  2316.       rw = NI - 1;        /* low guard word */
  2317.       rmsk = 0xffff;
  2318.       rmbit = 0x8000;
  2319.       re = rw - 1;
  2320.       rebit = 1;
  2321.       break;
  2322.     case 113:
  2323.       rw = 10;
  2324.       rmsk = 0x7fff;
  2325.       rmbit = 0x4000;
  2326.       rebit = 0x8000;
  2327.       re = rw;
  2328.       break;
  2329.     case 64:
  2330.       rw = 7;
  2331.       rmsk = 0xffff;
  2332.       rmbit = 0x8000;
  2333.       re = rw - 1;
  2334.       rebit = 1;
  2335.       break;
  2336.       /* For DEC or IBM arithmetic */
  2337.     case 56:
  2338.       rw = 6;
  2339.       rmsk = 0xff;
  2340.       rmbit = 0x80;
  2341.       rebit = 0x100;
  2342.       re = rw;
  2343.       break;
  2344.     case 53:
  2345.       rw = 6;
  2346.       rmsk = 0x7ff;
  2347.       rmbit = 0x0400;
  2348.       rebit = 0x800;
  2349.       re = rw;
  2350.       break;
  2351.     case 24:
  2352.       rw = 4;
  2353.       rmsk = 0xff;
  2354.       rmbit = 0x80;
  2355.       rebit = 0x100;
  2356.       re = rw;
  2357.       break;
  2358.     }
  2359.       rbit[re] = rebit;
  2360.       rlast = rndprc;
  2361.     }
  2362.  
  2363.   /* Shift down 1 temporarily if the data structure has an implied
  2364.      most significant bit and the number is denormal.  */
  2365.   if ((exp <= 0) && (rndprc != 64) && (rndprc != NBITS))
  2366.     {
  2367.       lost |= s[NI - 1] & 1;
  2368.       eshdn1 (s);
  2369.     }
  2370.   /* Clear out all bits below the rounding bit,
  2371.      remembering in r if any were nonzero.  */
  2372.   r = s[rw] & rmsk;
  2373.   if (rndprc < NBITS)
  2374.     {
  2375.       i = rw + 1;
  2376.       while (i < NI)
  2377.     {
  2378.       if (s[i])
  2379.         r |= 1;
  2380.       s[i] = 0;
  2381.       ++i;
  2382.     }
  2383.     }
  2384.   s[rw] &= ~rmsk;
  2385.   if ((r & rmbit) != 0)
  2386.     {
  2387.       if (r == rmbit)
  2388.     {
  2389.       if (lost == 0)
  2390.         {            /* round to even */
  2391.           if ((s[re] & rebit) == 0)
  2392.         goto mddone;
  2393.         }
  2394.       else
  2395.         {
  2396.           if (subflg != 0)
  2397.         goto mddone;
  2398.         }
  2399.     }
  2400.       eaddm (rbit, s);
  2401.     }
  2402.  mddone:
  2403.   if ((exp <= 0) && (rndprc != 64) && (rndprc != NBITS))
  2404.     {
  2405.       eshup1 (s);
  2406.     }
  2407.   if (s[2] != 0)
  2408.     {                /* overflow on roundoff */
  2409.       eshdn1 (s);
  2410.       exp += 1;
  2411.     }
  2412.  mdfin:
  2413.   s[NI - 1] = 0;
  2414.   if (exp >= 32767L)
  2415.     {
  2416. #ifndef INFINITY
  2417.     overf:
  2418. #endif
  2419. #ifdef INFINITY
  2420.       s[1] = 32767;
  2421.       for (i = 2; i < NI - 1; i++)
  2422.     s[i] = 0;
  2423.       if (extra_warnings)
  2424.     warning ("floating point overflow");
  2425. #else
  2426.       s[1] = 32766;
  2427.       s[2] = 0;
  2428.       for (i = M + 1; i < NI - 1; i++)
  2429.     s[i] = 0xffff;
  2430.       s[NI - 1] = 0;
  2431.       if ((rndprc < 64) || (rndprc == 113))
  2432.     {
  2433.       s[rw] &= ~rmsk;
  2434.       if (rndprc == 24)
  2435.         {
  2436.           s[5] = 0;
  2437.           s[6] = 0;
  2438.         }
  2439.     }
  2440. #endif
  2441.       return;
  2442.     }
  2443.   if (exp < 0)
  2444.     s[1] = 0;
  2445.   else
  2446.     s[1] = (unsigned EMUSHORT) exp;
  2447. }
  2448.  
  2449.  
  2450.  
  2451. /*  Subtract external format numbers.  */
  2452.  
  2453. static int subflg = 0;
  2454.  
  2455. static void 
  2456. esub (a, b, c)
  2457.      unsigned EMUSHORT *a, *b, *c;
  2458. {
  2459.  
  2460. #ifdef NANS
  2461.   if (eisnan (a))
  2462.     {
  2463.       emov (a, c);
  2464.       return;
  2465.     }
  2466.   if (eisnan (b))
  2467.     {
  2468.       emov (b, c);
  2469.       return;
  2470.     }
  2471. /* Infinity minus infinity is a NaN.
  2472.    Test for subtracting infinities of the same sign. */
  2473.   if (eisinf (a) && eisinf (b)
  2474.       && ((eisneg (a) ^ eisneg (b)) == 0))
  2475.     {
  2476.       mtherr ("esub", INVALID);
  2477.       enan (c, 0);
  2478.       return;
  2479.     }
  2480. #endif
  2481.   subflg = 1;
  2482.   eadd1 (a, b, c);
  2483. }
  2484.  
  2485.  
  2486. /* Add.  */
  2487.  
  2488. static void 
  2489. eadd (a, b, c)
  2490.      unsigned EMUSHORT *a, *b, *c;
  2491. {
  2492.  
  2493. #ifdef NANS
  2494. /* NaN plus anything is a NaN. */
  2495.   if (eisnan (a))
  2496.     {
  2497.       emov (a, c);
  2498.       return;
  2499.     }
  2500.   if (eisnan (b))
  2501.     {
  2502.       emov (b, c);
  2503.       return;
  2504.     }
  2505. /* Infinity minus infinity is a NaN.
  2506.    Test for adding infinities of opposite signs. */
  2507.   if (eisinf (a) && eisinf (b)
  2508.       && ((eisneg (a) ^ eisneg (b)) != 0))
  2509.     {
  2510.       mtherr ("esub", INVALID);
  2511.       enan (c, 0);
  2512.       return;
  2513.     }
  2514. #endif
  2515.   subflg = 0;
  2516.   eadd1 (a, b, c);
  2517. }
  2518.  
  2519. static void 
  2520. eadd1 (a, b, c)
  2521.      unsigned EMUSHORT *a, *b, *c;
  2522. {
  2523.   unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
  2524.   int i, lost, j, k;
  2525.   EMULONG lt, lta, ltb;
  2526.  
  2527. #ifdef INFINITY
  2528.   if (eisinf (a))
  2529.     {
  2530.       emov (a, c);
  2531.       if (subflg)
  2532.     eneg (c);
  2533.       return;
  2534.     }
  2535.   if (eisinf (b))
  2536.     {
  2537.       emov (b, c);
  2538.       return;
  2539.     }
  2540. #endif
  2541.   emovi (a, ai);
  2542.   emovi (b, bi);
  2543.   if (subflg)
  2544.     ai[0] = ~ai[0];
  2545.  
  2546.   /* compare exponents */
  2547.   lta = ai[E];
  2548.   ltb = bi[E];
  2549.   lt = lta - ltb;
  2550.   if (lt > 0L)
  2551.     {                /* put the larger number in bi */
  2552.       emovz (bi, ci);
  2553.       emovz (ai, bi);
  2554.       emovz (ci, ai);
  2555.       ltb = bi[E];
  2556.       lt = -lt;
  2557.     }
  2558.   lost = 0;
  2559.   if (lt != 0L)
  2560.     {
  2561.       if (lt < (EMULONG) (-NBITS - 1))
  2562.     goto done;        /* answer same as larger addend */
  2563.       k = (int) lt;
  2564.       lost = eshift (ai, k);    /* shift the smaller number down */
  2565.     }
  2566.   else
  2567.     {
  2568.       /* exponents were the same, so must compare significands */
  2569.       i = ecmpm (ai, bi);
  2570.       if (i == 0)
  2571.     {            /* the numbers are identical in magnitude */
  2572.       /* if different signs, result is zero */
  2573.       if (ai[0] != bi[0])
  2574.         {
  2575.           eclear (c);
  2576.           return;
  2577.         }
  2578.       /* if same sign, result is double */
  2579.       /* double denomalized tiny number */
  2580.       if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
  2581.         {
  2582.           eshup1 (bi);
  2583.           goto done;
  2584.         }
  2585.       /* add 1 to exponent unless both are zero! */
  2586.       for (j = 1; j < NI - 1; j++)
  2587.         {
  2588.           if (bi[j] != 0)
  2589.         {
  2590.           /* This could overflow, but let emovo take care of that. */
  2591.           ltb += 1;
  2592.           break;
  2593.         }
  2594.         }
  2595.       bi[E] = (unsigned EMUSHORT) ltb;
  2596.       goto done;
  2597.     }
  2598.       if (i > 0)
  2599.     {            /* put the larger number in bi */
  2600.       emovz (bi, ci);
  2601.       emovz (ai, bi);
  2602.       emovz (ci, ai);
  2603.     }
  2604.     }
  2605.   if (ai[0] == bi[0])
  2606.     {
  2607.       eaddm (ai, bi);
  2608.       subflg = 0;
  2609.     }
  2610.   else
  2611.     {
  2612.       esubm (ai, bi);
  2613.       subflg = 1;
  2614.     }
  2615.   emdnorm (bi, lost, subflg, ltb, 64);
  2616.  
  2617.  done:
  2618.   emovo (bi, c);
  2619. }
  2620.  
  2621.  
  2622.  
  2623. /* Divide.  */
  2624.  
  2625. static void 
  2626. ediv (a, b, c)
  2627.      unsigned EMUSHORT *a, *b, *c;
  2628. {
  2629.   unsigned EMUSHORT ai[NI], bi[NI];
  2630.   int i;
  2631.   EMULONG lt, lta, ltb;
  2632.  
  2633. #ifdef NANS
  2634. /* Return any NaN input. */
  2635.   if (eisnan (a))
  2636.     {
  2637.     emov (a, c);
  2638.     return;
  2639.     }
  2640.   if (eisnan (b))
  2641.     {
  2642.     emov (b, c);
  2643.     return;
  2644.     }
  2645. /* Zero over zero, or infinity over infinity, is a NaN. */
  2646.   if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
  2647.       || (eisinf (a) && eisinf (b)))
  2648.     {
  2649.     mtherr ("ediv", INVALID);
  2650.     enan (c, eisneg (a) ^ eisneg (b));
  2651.     return;
  2652.     }
  2653. #endif
  2654. /* Infinity over anything else is infinity. */
  2655. #ifdef INFINITY
  2656.   if (eisinf (b))
  2657.     {
  2658.       if (eisneg (a) ^ eisneg (b))
  2659.     *(c + (NE - 1)) = 0x8000;
  2660.       else
  2661.     *(c + (NE - 1)) = 0;
  2662.       einfin (c);
  2663.       return;
  2664.     }
  2665. /* Anything else over infinity is zero. */
  2666.   if (eisinf (a))
  2667.     {
  2668.       eclear (c);
  2669.       return;
  2670.     }
  2671. #endif
  2672.   emovi (a, ai);
  2673.   emovi (b, bi);
  2674.   lta = ai[E];
  2675.   ltb = bi[E];
  2676.   if (bi[E] == 0)
  2677.     {                /* See if numerator is zero. */
  2678.       for (i = 1; i < NI - 1; i++)
  2679.     {
  2680.       if (bi[i] != 0)
  2681.         {
  2682.           ltb -= enormlz (bi);
  2683.           goto dnzro1;
  2684.         }
  2685.     }
  2686.       eclear (c);
  2687.       return;
  2688.     }
  2689.  dnzro1:
  2690.  
  2691.   if (ai[E] == 0)
  2692.     {                /* possible divide by zero */
  2693.       for (i = 1; i < NI - 1; i++)
  2694.     {
  2695.       if (ai[i] != 0)
  2696.         {
  2697.           lta -= enormlz (ai);
  2698.           goto dnzro2;
  2699.         }
  2700.     }
  2701.       if (ai[0] == bi[0])
  2702.     *(c + (NE - 1)) = 0;
  2703.       else
  2704.     *(c + (NE - 1)) = 0x8000;
  2705. /* Divide by zero is not an invalid operation.
  2706.    It is a divide-by-zero operation!   */
  2707.       einfin (c);
  2708.       mtherr ("ediv", SING);
  2709.       return;
  2710.     }
  2711.  dnzro2:
  2712.  
  2713.   i = edivm (ai, bi);
  2714.   /* calculate exponent */
  2715.   lt = ltb - lta + EXONE;
  2716.   emdnorm (bi, i, 0, lt, 64);
  2717.   /* set the sign */
  2718.   if (ai[0] == bi[0])
  2719.     bi[0] = 0;
  2720.   else
  2721.     bi[0] = 0Xffff;
  2722.   emovo (bi, c);
  2723. }
  2724.  
  2725.  
  2726.  
  2727. /* Multiply.  */
  2728.  
  2729. static void 
  2730. emul (a, b, c)
  2731.      unsigned EMUSHORT *a, *b, *c;
  2732. {
  2733.   unsigned EMUSHORT ai[NI], bi[NI];
  2734.   int i, j;
  2735.   EMULONG lt, lta, ltb;
  2736.  
  2737. #ifdef NANS
  2738. /* NaN times anything is the same NaN. */
  2739.   if (eisnan (a))
  2740.     {
  2741.     emov (a, c);
  2742.     return;
  2743.     }
  2744.   if (eisnan (b))
  2745.     {
  2746.     emov (b, c);
  2747.     return;
  2748.     }
  2749. /* Zero times infinity is a NaN. */
  2750.   if ((eisinf (a) && (ecmp (b, ezero) == 0))
  2751.       || (eisinf (b) && (ecmp (a, ezero) == 0)))
  2752.     {
  2753.     mtherr ("emul", INVALID);
  2754.     enan (c, eisneg (a) ^ eisneg (b));
  2755.     return;
  2756.     }
  2757. #endif
  2758. /* Infinity times anything else is infinity. */
  2759. #ifdef INFINITY
  2760.   if (eisinf (a) || eisinf (b))
  2761.     {
  2762.       if (eisneg (a) ^ eisneg (b))
  2763.     *(c + (NE - 1)) = 0x8000;
  2764.       else
  2765.     *(c + (NE - 1)) = 0;
  2766.       einfin (c);
  2767.       return;
  2768.     }
  2769. #endif
  2770.   emovi (a, ai);
  2771.   emovi (b, bi);
  2772.   lta = ai[E];
  2773.   ltb = bi[E];
  2774.   if (ai[E] == 0)
  2775.     {
  2776.       for (i = 1; i < NI - 1; i++)
  2777.     {
  2778.       if (ai[i] != 0)
  2779.         {
  2780.           lta -= enormlz (ai);
  2781.           goto mnzer1;
  2782.         }
  2783.     }
  2784.       eclear (c);
  2785.       return;
  2786.     }
  2787.  mnzer1:
  2788.  
  2789.   if (bi[E] == 0)
  2790.     {
  2791.       for (i = 1; i < NI - 1; i++)
  2792.     {
  2793.       if (bi[i] != 0)
  2794.         {
  2795.           ltb -= enormlz (bi);
  2796.           goto mnzer2;
  2797.         }
  2798.     }
  2799.       eclear (c);
  2800.       return;
  2801.     }
  2802.  mnzer2:
  2803.  
  2804.   /* Multiply significands */
  2805.   j = emulm (ai, bi);
  2806.   /* calculate exponent */
  2807.   lt = lta + ltb - (EXONE - 1);
  2808.   emdnorm (bi, j, 0, lt, 64);
  2809.   /* calculate sign of product */
  2810.   if (ai[0] == bi[0])
  2811.     bi[0] = 0;
  2812.   else
  2813.     bi[0] = 0xffff;
  2814.   emovo (bi, c);
  2815. }
  2816.  
  2817.  
  2818.  
  2819.  
  2820. /* Convert IEEE double precision to e type.  */
  2821.  
  2822. static void
  2823. e53toe (pe, y)
  2824.      unsigned EMUSHORT *pe, *y;
  2825. {
  2826. #ifdef DEC
  2827.  
  2828.   dectoe (pe, y);        /* see etodec.c */
  2829.  
  2830. #else
  2831. #ifdef IBM
  2832.  
  2833.   ibmtoe (pe, y, DFmode);
  2834.  
  2835. #else
  2836.   register unsigned EMUSHORT r;
  2837.   register unsigned EMUSHORT *e, *p;
  2838.   unsigned EMUSHORT yy[NI];
  2839.   int denorm, k;
  2840.  
  2841.   e = pe;
  2842.   denorm = 0;            /* flag if denormalized number */
  2843.   ecleaz (yy);
  2844. #ifdef IBMPC
  2845.   e += 3;
  2846. #endif
  2847.   r = *e;
  2848.   yy[0] = 0;
  2849.   if (r & 0x8000)
  2850.     yy[0] = 0xffff;
  2851.   yy[M] = (r & 0x0f) | 0x10;
  2852.   r &= ~0x800f;            /* strip sign and 4 significand bits */
  2853. #ifdef INFINITY
  2854.   if (r == 0x7ff0)
  2855.     {
  2856. #ifdef NANS
  2857. #ifdef IBMPC
  2858.       if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
  2859.       || (pe[1] != 0) || (pe[0] != 0))
  2860.     {
  2861.       enan (y, yy[0] != 0);
  2862.       return;
  2863.     }
  2864. #else
  2865.       if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
  2866.       || (pe[2] != 0) || (pe[3] != 0))
  2867.     {
  2868.       enan (y, yy[0] != 0);
  2869.       return;
  2870.     }
  2871. #endif
  2872. #endif  /* NANS */
  2873.       eclear (y);
  2874.       einfin (y);
  2875.       if (yy[0])
  2876.     eneg (y);
  2877.       return;
  2878.     }
  2879. #endif  /* INFINITY */
  2880.   r >>= 4;
  2881.   /* If zero exponent, then the significand is denormalized.
  2882.      So take back the understood high significand bit. */
  2883.  
  2884.   if (r == 0)
  2885.     {
  2886.       denorm = 1;
  2887.       yy[M] &= ~0x10;
  2888.     }
  2889.   r += EXONE - 01777;
  2890.   yy[E] = r;
  2891.   p = &yy[M + 1];
  2892. #ifdef IBMPC
  2893.   *p++ = *(--e);
  2894.   *p++ = *(--e);
  2895.   *p++ = *(--e);
  2896. #endif
  2897. #ifdef MIEEE
  2898.   ++e;
  2899.   *p++ = *e++;
  2900.   *p++ = *e++;
  2901.   *p++ = *e++;
  2902. #endif
  2903.   eshift (yy, -5);
  2904.   if (denorm)
  2905.     {                /* if zero exponent, then normalize the significand */
  2906.       if ((k = enormlz (yy)) > NBITS)
  2907.     ecleazs (yy);
  2908.       else
  2909.     yy[E] -= (unsigned EMUSHORT) (k - 1);
  2910.     }
  2911.   emovo (yy, y);
  2912. #endif /* not IBM */
  2913. #endif /* not DEC */
  2914. }
  2915.  
  2916. static void 
  2917. e64toe (pe, y)
  2918.      unsigned EMUSHORT *pe, *y;
  2919. {
  2920.   unsigned EMUSHORT yy[NI];
  2921.   unsigned EMUSHORT *e, *p, *q;
  2922.   int i;
  2923.  
  2924.   e = pe;
  2925.   p = yy;
  2926.   for (i = 0; i < NE - 5; i++)
  2927.     *p++ = 0;
  2928. #ifdef IBMPC
  2929.   for (i = 0; i < 5; i++)
  2930.     *p++ = *e++;
  2931. #endif
  2932. /* This precision is not ordinarily supported on DEC or IBM. */
  2933. #ifdef DEC
  2934.   for (i = 0; i < 5; i++)
  2935.     *p++ = *e++;
  2936. #endif
  2937. #ifdef IBM
  2938.   p = &yy[0] + (NE - 1);
  2939.   *p-- = *e++;
  2940.   ++e;
  2941.   for (i = 0; i < 5; i++)
  2942.     *p-- = *e++;
  2943. #endif
  2944. #ifdef MIEEE
  2945.   p = &yy[0] + (NE - 1);
  2946.   *p-- = *e++;
  2947.   ++e;
  2948.   for (i = 0; i < 4; i++)
  2949.     *p-- = *e++;
  2950. #endif
  2951.   p = yy;
  2952.   q = y;
  2953. #ifdef INFINITY
  2954.   if (*p == 0x7fff)
  2955.     {
  2956. #ifdef NANS
  2957. #ifdef IBMPC
  2958.       for (i = 0; i < 4; i++)
  2959.     {
  2960.       if (pe[i] != 0)
  2961.         {
  2962.           enan (y, (*p & 0x8000) != 0);
  2963.           return;
  2964.         }
  2965.     }
  2966. #else
  2967.       for (i = 1; i <= 4; i++)
  2968.     {
  2969.       if (pe[i] != 0)
  2970.         {
  2971.           enan (y, (*p & 0x8000) != 0);
  2972.           return;
  2973.         }
  2974.     }
  2975. #endif
  2976. #endif /* NANS */
  2977.       eclear (y);
  2978.       einfin (y);
  2979.       if (*p & 0x8000)
  2980.     eneg (y);
  2981.       return;
  2982.     }
  2983. #endif  /* INFINITY */
  2984.   for (i = 0; i < NE; i++)
  2985.     *q++ = *p++;
  2986. }
  2987.  
  2988.  
  2989. static void 
  2990. e113toe (pe, y)
  2991.      unsigned EMUSHORT *pe, *y;
  2992. {
  2993.   register unsigned EMUSHORT r;
  2994.   unsigned EMUSHORT *e, *p;
  2995.   unsigned EMUSHORT yy[NI];
  2996.   int denorm, i;
  2997.  
  2998.   e = pe;
  2999.   denorm = 0;
  3000.   ecleaz (yy);
  3001. #ifdef IBMPC
  3002.   e += 7;
  3003. #endif
  3004.   r = *e;
  3005.   yy[0] = 0;
  3006.   if (r & 0x8000)
  3007.     yy[0] = 0xffff;
  3008.   r &= 0x7fff;
  3009. #ifdef INFINITY
  3010.   if (r == 0x7fff)
  3011.     {
  3012. #ifdef NANS
  3013. #ifdef IBMPC
  3014.       for (i = 0; i < 7; i++)
  3015.     {
  3016.       if (pe[i] != 0)
  3017.         {
  3018.           enan (y, yy[0] != 0);
  3019.           return;
  3020.         }
  3021.     }
  3022. #else
  3023.       for (i = 1; i < 8; i++)
  3024.     {
  3025.       if (pe[i] != 0)
  3026.         {
  3027.           enan (y, yy[0] != 0);
  3028.           return;
  3029.         }
  3030.     }
  3031. #endif
  3032. #endif /* NANS */
  3033.       eclear (y);
  3034.       einfin (y);
  3035.       if (yy[0])
  3036.     eneg (y);
  3037.       return;
  3038.     }
  3039. #endif  /* INFINITY */
  3040.   yy[E] = r;
  3041.   p = &yy[M + 1];
  3042. #ifdef IBMPC
  3043.   for (i = 0; i < 7; i++)
  3044.     *p++ = *(--e);
  3045. #endif
  3046. #ifdef MIEEE
  3047.   ++e;
  3048.   for (i = 0; i < 7; i++)
  3049.     *p++ = *e++;
  3050. #endif
  3051. /* If denormal, remove the implied bit; else shift down 1. */
  3052.   if (r == 0)
  3053.     {
  3054.       yy[M] = 0;
  3055.     }
  3056.   else
  3057.     {
  3058.       yy[M] = 1;
  3059.       eshift (yy, -1);
  3060.     }
  3061.   emovo (yy, y);
  3062. }
  3063.  
  3064.  
  3065. /* Convert IEEE single precision to e type.  */
  3066.  
  3067. static void 
  3068. e24toe (pe, y)
  3069.      unsigned EMUSHORT *pe, *y;
  3070. {
  3071. #ifdef IBM
  3072.  
  3073.   ibmtoe (pe, y, SFmode);
  3074.  
  3075. #else
  3076.   register unsigned EMUSHORT r;
  3077.   register unsigned EMUSHORT *e, *p;
  3078.   unsigned EMUSHORT yy[NI];
  3079.   int denorm, k;
  3080.  
  3081.   e = pe;
  3082.   denorm = 0;            /* flag if denormalized number */
  3083.   ecleaz (yy);
  3084. #ifdef IBMPC
  3085.   e += 1;
  3086. #endif
  3087. #ifdef DEC
  3088.   e += 1;
  3089. #endif
  3090.   r = *e;
  3091.   yy[0] = 0;
  3092.   if (r & 0x8000)
  3093.     yy[0] = 0xffff;
  3094.   yy[M] = (r & 0x7f) | 0200;
  3095.   r &= ~0x807f;            /* strip sign and 7 significand bits */
  3096. #ifdef INFINITY
  3097.   if (r == 0x7f80)
  3098.     {
  3099. #ifdef NANS
  3100. #ifdef MIEEE
  3101.       if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
  3102.     {
  3103.       enan (y, yy[0] != 0);
  3104.       return;
  3105.     }
  3106. #else
  3107.       if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
  3108.     {
  3109.       enan (y, yy[0] != 0);
  3110.       return;
  3111.     }
  3112. #endif
  3113. #endif  /* NANS */
  3114.       eclear (y);
  3115.       einfin (y);
  3116.       if (yy[0])
  3117.     eneg (y);
  3118.       return;
  3119.     }
  3120. #endif  /* INFINITY */
  3121.   r >>= 7;
  3122.   /* If zero exponent, then the significand is denormalized.
  3123.      So take back the understood high significand bit. */
  3124.   if (r == 0)
  3125.     {
  3126.       denorm = 1;
  3127.       yy[M] &= ~0200;
  3128.     }
  3129.   r += EXONE - 0177;
  3130.   yy[E] = r;
  3131.   p = &yy[M + 1];
  3132. #ifdef IBMPC
  3133.   *p++ = *(--e);
  3134. #endif
  3135. #ifdef DEC
  3136.   *p++ = *(--e);
  3137. #endif
  3138. #ifdef MIEEE
  3139.   ++e;
  3140.   *p++ = *e++;
  3141. #endif
  3142.   eshift (yy, -8);
  3143.   if (denorm)
  3144.     {                /* if zero exponent, then normalize the significand */
  3145.       if ((k = enormlz (yy)) > NBITS)
  3146.     ecleazs (yy);
  3147.       else
  3148.     yy[E] -= (unsigned EMUSHORT) (k - 1);
  3149.     }
  3150.   emovo (yy, y);
  3151. #endif /* not IBM */
  3152. }
  3153.  
  3154.  
  3155. static void 
  3156. etoe113 (x, e)
  3157.      unsigned EMUSHORT *x, *e;
  3158. {
  3159.   unsigned EMUSHORT xi[NI];
  3160.   EMULONG exp;
  3161.   int rndsav;
  3162.  
  3163. #ifdef NANS
  3164.   if (eisnan (x))
  3165.     {
  3166.       make_nan (e, eisneg (x), TFmode);
  3167.       return;
  3168.     }
  3169. #endif
  3170.   emovi (x, xi);
  3171.   exp = (EMULONG) xi[E];
  3172. #ifdef INFINITY
  3173.   if (eisinf (x))
  3174.     goto nonorm;
  3175. #endif
  3176.   /* round off to nearest or even */
  3177.   rndsav = rndprc;
  3178.   rndprc = 113;
  3179.   emdnorm (xi, 0, 0, exp, 64);
  3180.   rndprc = rndsav;
  3181.  nonorm:
  3182.   toe113 (xi, e);
  3183. }
  3184.  
  3185. /* Move out internal format to ieee long double */
  3186.  
  3187. static void 
  3188. toe113 (a, b)
  3189.      unsigned EMUSHORT *a, *b;
  3190. {
  3191.   register unsigned EMUSHORT *p, *q;
  3192.   unsigned EMUSHORT i;
  3193.  
  3194. #ifdef NANS
  3195.   if (eiisnan (a))
  3196.     {
  3197.       make_nan (b, eiisneg (a), TFmode);
  3198.       return;
  3199.     }
  3200. #endif
  3201.   p = a;
  3202. #ifdef MIEEE
  3203.   q = b;
  3204. #else
  3205.   q = b + 7;            /* point to output exponent */
  3206. #endif
  3207.  
  3208.   /* If not denormal, delete the implied bit. */
  3209.   if (a[E] != 0)
  3210.     {
  3211.       eshup1 (a);
  3212.     }
  3213.   /* combine sign and exponent */
  3214.   i = *p++;
  3215. #ifdef MIEEE
  3216.   if (i)
  3217.     *q++ = *p++ | 0x8000;
  3218.   else
  3219.     *q++ = *p++;
  3220. #else
  3221.   if (i)
  3222.     *q-- = *p++ | 0x8000;
  3223.   else
  3224.     *q-- = *p++;
  3225. #endif
  3226.   /* skip over guard word */
  3227.   ++p;
  3228.   /* move the significand */
  3229. #ifdef MIEEE
  3230.   for (i = 0; i < 7; i++)
  3231.     *q++ = *p++;
  3232. #else
  3233.   for (i = 0; i < 7; i++)
  3234.     *q-- = *p++;
  3235. #endif
  3236. }
  3237.  
  3238. static void 
  3239. etoe64 (x, e)
  3240.      unsigned EMUSHORT *x, *e;
  3241. {
  3242.   unsigned EMUSHORT xi[NI];
  3243.   EMULONG exp;
  3244.   int rndsav;
  3245.  
  3246. #ifdef NANS
  3247.   if (eisnan (x))
  3248.     {
  3249.       make_nan (e, eisneg (x), XFmode);
  3250.       return;
  3251.     }
  3252. #endif
  3253.   emovi (x, xi);
  3254.   /* adjust exponent for offset */
  3255.   exp = (EMULONG) xi[E];
  3256. #ifdef INFINITY
  3257.   if (eisinf (x))
  3258.     goto nonorm;
  3259. #endif
  3260.   /* round off to nearest or even */
  3261.   rndsav = rndprc;
  3262.   rndprc = 64;
  3263.   emdnorm (xi, 0, 0, exp, 64);
  3264.   rndprc = rndsav;
  3265.  nonorm:
  3266.   toe64 (xi, e);
  3267. }
  3268.  
  3269.  
  3270. /* Move out internal format to ieee long double. */
  3271.  
  3272. static void 
  3273. toe64 (a, b)
  3274.      unsigned EMUSHORT *a, *b;
  3275. {
  3276.   register unsigned EMUSHORT *p, *q;
  3277.   unsigned EMUSHORT i;
  3278.  
  3279. #ifdef NANS
  3280.   if (eiisnan (a))
  3281.     {
  3282.       make_nan (b, eiisneg (a), XFmode);
  3283.       return;
  3284.     }
  3285. #endif
  3286.   p = a;
  3287. #if defined(MIEEE) || defined(IBM)
  3288.   q = b;
  3289. #else
  3290.   q = b + 4;            /* point to output exponent */
  3291. #if LONG_DOUBLE_TYPE_SIZE == 96
  3292.   /* Clear the last two bytes of 12-byte Intel format */
  3293.   *(q+1) = 0;
  3294. #endif
  3295. #endif
  3296.  
  3297.   /* combine sign and exponent */
  3298.   i = *p++;
  3299. #if defined(MIEEE) || defined(IBM)
  3300.   if (i)
  3301.     *q++ = *p++ | 0x8000;
  3302.   else
  3303.     *q++ = *p++;
  3304.   *q++ = 0;
  3305. #else
  3306.   if (i)
  3307.     *q-- = *p++ | 0x8000;
  3308.   else
  3309.     *q-- = *p++;
  3310. #endif
  3311.   /* skip over guard word */
  3312.   ++p;
  3313.   /* move the significand */
  3314. #if defined(MIEEE) || defined(IBM)
  3315.   for (i = 0; i < 4; i++)
  3316.     *q++ = *p++;
  3317. #else
  3318.   for (i = 0; i < 4; i++)
  3319.     *q-- = *p++;
  3320. #endif
  3321. }
  3322.  
  3323.  
  3324. /* e type to IEEE double precision.  */
  3325.  
  3326. #ifdef DEC
  3327.  
  3328. static void 
  3329. etoe53 (x, e)
  3330.      unsigned EMUSHORT *x, *e;
  3331. {
  3332.   etodec (x, e);        /* see etodec.c */
  3333. }
  3334.  
  3335. static void 
  3336. toe53 (x, y)
  3337.      unsigned EMUSHORT *x, *y;
  3338. {
  3339.   todec (x, y);
  3340. }
  3341.  
  3342. #else
  3343. #ifdef IBM
  3344.  
  3345. static void 
  3346. etoe53 (x, e)
  3347.      unsigned EMUSHORT *x, *e;
  3348. {
  3349.   etoibm (x, e, DFmode);
  3350. }
  3351.  
  3352. static void 
  3353. toe53 (x, y)
  3354.      unsigned EMUSHORT *x, *y;
  3355. {
  3356.   toibm (x, y, DFmode);
  3357. }
  3358.  
  3359. #else  /* it's neither DEC nor IBM */
  3360.  
  3361. static void 
  3362. etoe53 (x, e)
  3363.      unsigned EMUSHORT *x, *e;
  3364. {
  3365.   unsigned EMUSHORT xi[NI];
  3366.   EMULONG exp;
  3367.   int rndsav;
  3368.  
  3369. #ifdef NANS
  3370.   if (eisnan (x))
  3371.     {
  3372.       make_nan (e, eisneg (x), DFmode);
  3373.       return;
  3374.     }
  3375. #endif
  3376.   emovi (x, xi);
  3377.   /* adjust exponent for offsets */
  3378.   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
  3379. #ifdef INFINITY
  3380.   if (eisinf (x))
  3381.     goto nonorm;
  3382. #endif
  3383.   /* round off to nearest or even */
  3384.   rndsav = rndprc;
  3385.   rndprc = 53;
  3386.   emdnorm (xi, 0, 0, exp, 64);
  3387.   rndprc = rndsav;
  3388.  nonorm:
  3389.   toe53 (xi, e);
  3390. }
  3391.  
  3392.  
  3393. static void 
  3394. toe53 (x, y)
  3395.      unsigned EMUSHORT *x, *y;
  3396. {
  3397.   unsigned EMUSHORT i;
  3398.   unsigned EMUSHORT *p;
  3399.  
  3400. #ifdef NANS
  3401.   if (eiisnan (x))
  3402.     {
  3403.       make_nan (y, eiisneg (x), DFmode);
  3404.       return;
  3405.     }
  3406. #endif
  3407.   p = &x[0];
  3408. #ifdef IBMPC
  3409.   y += 3;
  3410. #endif
  3411.   *y = 0;            /* output high order */
  3412.   if (*p++)
  3413.     *y = 0x8000;        /* output sign bit */
  3414.  
  3415.   i = *p++;
  3416.   if (i >= (unsigned int) 2047)
  3417.     {                /* Saturate at largest number less than infinity. */
  3418. #ifdef INFINITY
  3419.       *y |= 0x7ff0;
  3420. #ifdef IBMPC
  3421.       *(--y) = 0;
  3422.       *(--y) = 0;
  3423.       *(--y) = 0;
  3424. #endif
  3425. #ifdef MIEEE
  3426.       ++y;
  3427.       *y++ = 0;
  3428.       *y++ = 0;
  3429.       *y++ = 0;
  3430. #endif
  3431. #else
  3432.       *y |= (unsigned EMUSHORT) 0x7fef;
  3433. #ifdef IBMPC
  3434.       *(--y) = 0xffff;
  3435.       *(--y) = 0xffff;
  3436.       *(--y) = 0xffff;
  3437. #endif
  3438. #ifdef MIEEE
  3439.       ++y;
  3440.       *y++ = 0xffff;
  3441.       *y++ = 0xffff;
  3442.       *y++ = 0xffff;
  3443. #endif
  3444. #endif
  3445.       return;
  3446.     }
  3447.   if (i == 0)
  3448.     {
  3449.       eshift (x, 4);
  3450.     }
  3451.   else
  3452.     {
  3453.       i <<= 4;
  3454.       eshift (x, 5);
  3455.     }
  3456.   i |= *p++ & (unsigned EMUSHORT) 0x0f;    /* *p = xi[M] */
  3457.   *y |= (unsigned EMUSHORT) i;    /* high order output already has sign bit set */
  3458. #ifdef IBMPC
  3459.   *(--y) = *p++;
  3460.   *(--y) = *p++;
  3461.   *(--y) = *p;
  3462. #endif
  3463. #ifdef MIEEE
  3464.   ++y;
  3465.   *y++ = *p++;
  3466.   *y++ = *p++;
  3467.   *y++ = *p++;
  3468. #endif
  3469. }
  3470.  
  3471. #endif /* not IBM */
  3472. #endif /* not DEC */
  3473.  
  3474.  
  3475.  
  3476. /* e type to IEEE single precision.  */
  3477.  
  3478. #ifdef IBM
  3479.  
  3480. static void 
  3481. etoe24 (x, e)
  3482.      unsigned EMUSHORT *x, *e;
  3483. {
  3484.   etoibm (x, e, SFmode);
  3485. }
  3486.  
  3487. static void 
  3488. toe24 (x, y)
  3489.      unsigned EMUSHORT *x, *y;
  3490. {
  3491.   toibm (x, y, SFmode);
  3492. }
  3493.  
  3494. #else
  3495.  
  3496. static void 
  3497. etoe24 (x, e)
  3498.      unsigned EMUSHORT *x, *e;
  3499. {
  3500.   EMULONG exp;
  3501.   unsigned EMUSHORT xi[NI];
  3502.   int rndsav;
  3503.  
  3504. #ifdef NANS
  3505.   if (eisnan (x))
  3506.     {
  3507.       make_nan (e, eisneg (x), SFmode);
  3508.       return;
  3509.     }
  3510. #endif
  3511.   emovi (x, xi);
  3512.   /* adjust exponent for offsets */
  3513.   exp = (EMULONG) xi[E] - (EXONE - 0177);
  3514. #ifdef INFINITY
  3515.   if (eisinf (x))
  3516.     goto nonorm;
  3517. #endif
  3518.   /* round off to nearest or even */
  3519.   rndsav = rndprc;
  3520.   rndprc = 24;
  3521.   emdnorm (xi, 0, 0, exp, 64);
  3522.   rndprc = rndsav;
  3523.  nonorm:
  3524.   toe24 (xi, e);
  3525. }
  3526.  
  3527. static void 
  3528. toe24 (x, y)
  3529.      unsigned EMUSHORT *x, *y;
  3530. {
  3531.   unsigned EMUSHORT i;
  3532.   unsigned EMUSHORT *p;
  3533.  
  3534. #ifdef NANS
  3535.   if (eiisnan (x))
  3536.     {
  3537.       make_nan (y, eiisneg (x), SFmode);
  3538.       return;
  3539.     }
  3540. #endif
  3541.   p = &x[0];
  3542. #ifdef IBMPC
  3543.   y += 1;
  3544. #endif
  3545. #ifdef DEC
  3546.   y += 1;
  3547. #endif
  3548.   *y = 0;            /* output high order */
  3549.   if (*p++)
  3550.     *y = 0x8000;        /* output sign bit */
  3551.  
  3552.   i = *p++;
  3553. /* Handle overflow cases. */
  3554.   if (i >= 255)
  3555.     {
  3556. #ifdef INFINITY
  3557.       *y |= (unsigned EMUSHORT) 0x7f80;
  3558. #ifdef IBMPC
  3559.       *(--y) = 0;
  3560. #endif
  3561. #ifdef DEC
  3562.       *(--y) = 0;
  3563. #endif
  3564. #ifdef MIEEE
  3565.       ++y;
  3566.       *y = 0;
  3567. #endif
  3568. #else  /* no INFINITY */
  3569.       *y |= (unsigned EMUSHORT) 0x7f7f;
  3570. #ifdef IBMPC
  3571.       *(--y) = 0xffff;
  3572. #endif
  3573. #ifdef DEC
  3574.       *(--y) = 0xffff;
  3575. #endif
  3576. #ifdef MIEEE
  3577.       ++y;
  3578.       *y = 0xffff;
  3579. #endif
  3580. #ifdef ERANGE
  3581.       errno = ERANGE;
  3582. #endif
  3583. #endif  /* no INFINITY */
  3584.       return;
  3585.     }
  3586.   if (i == 0)
  3587.     {
  3588.       eshift (x, 7);
  3589.     }
  3590.   else
  3591.     {
  3592.       i <<= 7;
  3593.       eshift (x, 8);
  3594.     }
  3595.   i |= *p++ & (unsigned EMUSHORT) 0x7f;    /* *p = xi[M] */
  3596.   *y |= i;            /* high order output already has sign bit set */
  3597. #ifdef IBMPC
  3598.   *(--y) = *p;
  3599. #endif
  3600. #ifdef DEC
  3601.   *(--y) = *p;
  3602. #endif
  3603. #ifdef MIEEE
  3604.   ++y;
  3605.   *y = *p;
  3606. #endif
  3607. }
  3608. #endif  /* not IBM */
  3609.  
  3610. /* Compare two e type numbers. 
  3611.    Return +1 if a > b
  3612.            0 if a == b
  3613.           -1 if a < b
  3614.           -2 if either a or b is a NaN.  */
  3615.  
  3616. static int 
  3617. ecmp (a, b)
  3618.      unsigned EMUSHORT *a, *b;
  3619. {
  3620.   unsigned EMUSHORT ai[NI], bi[NI];
  3621.   register unsigned EMUSHORT *p, *q;
  3622.   register int i;
  3623.   int msign;
  3624.  
  3625. #ifdef NANS
  3626.   if (eisnan (a)  || eisnan (b))
  3627.       return (-2);
  3628. #endif
  3629.   emovi (a, ai);
  3630.   p = ai;
  3631.   emovi (b, bi);
  3632.   q = bi;
  3633.  
  3634.   if (*p != *q)
  3635.     {                /* the signs are different */
  3636.       /* -0 equals + 0 */
  3637.       for (i = 1; i < NI - 1; i++)
  3638.     {
  3639.       if (ai[i] != 0)
  3640.         goto nzro;
  3641.       if (bi[i] != 0)
  3642.         goto nzro;
  3643.     }
  3644.       return (0);
  3645.     nzro:
  3646.       if (*p == 0)
  3647.     return (1);
  3648.       else
  3649.     return (-1);
  3650.     }
  3651.   /* both are the same sign */
  3652.   if (*p == 0)
  3653.     msign = 1;
  3654.   else
  3655.     msign = -1;
  3656.   i = NI - 1;
  3657.   do
  3658.     {
  3659.       if (*p++ != *q++)
  3660.     {
  3661.       goto diff;
  3662.     }
  3663.     }
  3664.   while (--i > 0);
  3665.  
  3666.   return (0);            /* equality */
  3667.  
  3668.  
  3669.  
  3670.  diff:
  3671.  
  3672.   if (*(--p) > *(--q))
  3673.     return (msign);        /* p is bigger */
  3674.   else
  3675.     return (-msign);        /* p is littler */
  3676. }
  3677.  
  3678.  
  3679.  
  3680.  
  3681. /* Find nearest integer to x = floor (x + 0.5).  */
  3682.  
  3683. static void 
  3684. eround (x, y)
  3685.      unsigned EMUSHORT *x, *y;
  3686. {
  3687.   eadd (ehalf, x, y);
  3688.   efloor (y, y);
  3689. }
  3690.  
  3691.  
  3692.  
  3693.  
  3694. /* Convert HOST_WIDE_INT to e type.  */
  3695.  
  3696. static void 
  3697. ltoe (lp, y)
  3698.      HOST_WIDE_INT *lp;
  3699.      unsigned EMUSHORT *y;
  3700. {
  3701.   unsigned EMUSHORT yi[NI];
  3702.   unsigned HOST_WIDE_INT ll;
  3703.   int k;
  3704.  
  3705.   ecleaz (yi);
  3706.   if (*lp < 0)
  3707.     {
  3708.       /* make it positive */
  3709.       ll = (unsigned HOST_WIDE_INT) (-(*lp));
  3710.       yi[0] = 0xffff;        /* put correct sign in the e type number */
  3711.     }
  3712.   else
  3713.     {
  3714.       ll = (unsigned HOST_WIDE_INT) (*lp);
  3715.     }
  3716.   /* move the long integer to yi significand area */
  3717. #if HOST_BITS_PER_WIDE_INT == 64
  3718.   yi[M] = (unsigned EMUSHORT) (ll >> 48);
  3719.   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
  3720.   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
  3721.   yi[M + 3] = (unsigned EMUSHORT) ll;
  3722.   yi[E] = EXONE + 47;        /* exponent if normalize shift count were 0 */
  3723. #else
  3724.   yi[M] = (unsigned EMUSHORT) (ll >> 16);
  3725.   yi[M + 1] = (unsigned EMUSHORT) ll;
  3726.   yi[E] = EXONE + 15;        /* exponent if normalize shift count were 0 */
  3727. #endif
  3728.  
  3729.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  3730.     ecleaz (yi);        /* it was zero */
  3731.   else
  3732.     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
  3733.   emovo (yi, y);        /* output the answer */
  3734. }
  3735.  
  3736. /* Convert unsigned HOST_WIDE_INT to e type.  */
  3737.  
  3738. static void 
  3739. ultoe (lp, y)
  3740.      unsigned HOST_WIDE_INT *lp;
  3741.      unsigned EMUSHORT *y;
  3742. {
  3743.   unsigned EMUSHORT yi[NI];
  3744.   unsigned HOST_WIDE_INT ll;
  3745.   int k;
  3746.  
  3747.   ecleaz (yi);
  3748.   ll = *lp;
  3749.  
  3750.   /* move the long integer to ayi significand area */
  3751. #if HOST_BITS_PER_WIDE_INT == 64
  3752.   yi[M] = (unsigned EMUSHORT) (ll >> 48);
  3753.   yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
  3754.   yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
  3755.   yi[M + 3] = (unsigned EMUSHORT) ll;
  3756.   yi[E] = EXONE + 47;        /* exponent if normalize shift count were 0 */
  3757. #else
  3758.   yi[M] = (unsigned EMUSHORT) (ll >> 16);
  3759.   yi[M + 1] = (unsigned EMUSHORT) ll;
  3760.   yi[E] = EXONE + 15;        /* exponent if normalize shift count were 0 */
  3761. #endif
  3762.  
  3763.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  3764.     ecleaz (yi);        /* it was zero */
  3765.   else
  3766.     yi[E] -= (unsigned EMUSHORT) k;  /* subtract shift count from exponent */
  3767.   emovo (yi, y);        /* output the answer */
  3768. }
  3769.  
  3770.  
  3771. /* Find signed HOST_WIDE_INT integer and floating point fractional
  3772.    parts of e-type (packed internal format) floating point input X.
  3773.    The integer output I has the sign of the input, except that
  3774.    positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
  3775.    The output e-type fraction FRAC is the positive fractional
  3776.    part of abs (X).  */
  3777.  
  3778. static void 
  3779. eifrac (x, i, frac)
  3780.      unsigned EMUSHORT *x;
  3781.      HOST_WIDE_INT *i;
  3782.      unsigned EMUSHORT *frac;
  3783. {
  3784.   unsigned EMUSHORT xi[NI];
  3785.   int j, k;
  3786.   unsigned HOST_WIDE_INT ll;
  3787.  
  3788.   emovi (x, xi);
  3789.   k = (int) xi[E] - (EXONE - 1);
  3790.   if (k <= 0)
  3791.     {
  3792.       /* if exponent <= 0, integer = 0 and real output is fraction */
  3793.       *i = 0L;
  3794.       emovo (xi, frac);
  3795.       return;
  3796.     }
  3797.   if (k > (HOST_BITS_PER_WIDE_INT - 1))
  3798.     {
  3799.       /* long integer overflow: output large integer
  3800.      and correct fraction  */
  3801.       if (xi[0])
  3802.     *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
  3803.       else
  3804.     {
  3805. #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
  3806.       /* In this case, let it overflow and convert as if unsigned.  */
  3807.       euifrac (x, &ll, frac);
  3808.       *i = (HOST_WIDE_INT) ll;
  3809.       return;
  3810. #else
  3811.       /* In other cases, return the largest positive integer.  */
  3812.       *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
  3813. #endif
  3814.     }
  3815.       eshift (xi, k);
  3816.       if (extra_warnings)
  3817.     warning ("overflow on truncation to integer");
  3818.     }
  3819.   else if (k > 16)
  3820.     {
  3821.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  3822.      then shift up by 16's.  */
  3823.       j = k - ((k >> 4) << 4);
  3824.       eshift (xi, j);
  3825.       ll = xi[M];
  3826.       k -= j;
  3827.       do
  3828.     {
  3829.       eshup6 (xi);
  3830.       ll = (ll << 16) | xi[M];
  3831.     }
  3832.       while ((k -= 16) > 0);
  3833.       *i = ll;
  3834.       if (xi[0])
  3835.     *i = -(*i);
  3836.     }
  3837.   else
  3838.       {
  3839.         /* shift not more than 16 bits */
  3840.           eshift (xi, k);
  3841.         *i = (HOST_WIDE_INT) xi[M] & 0xffff;
  3842.         if (xi[0])
  3843.       *i = -(*i);
  3844.       }
  3845.   xi[0] = 0;
  3846.   xi[E] = EXONE - 1;
  3847.   xi[M] = 0;
  3848.   if ((k = enormlz (xi)) > NBITS)
  3849.     ecleaz (xi);
  3850.   else
  3851.     xi[E] -= (unsigned EMUSHORT) k;
  3852.  
  3853.   emovo (xi, frac);
  3854. }
  3855.  
  3856.  
  3857. /* Find unsigned HOST_WIDE_INT integer and floating point fractional parts.
  3858.    A negative e type input yields integer output = 0
  3859.    but correct fraction.  */
  3860.  
  3861. static void 
  3862. euifrac (x, i, frac)
  3863.      unsigned EMUSHORT *x;
  3864.      unsigned HOST_WIDE_INT *i;
  3865.      unsigned EMUSHORT *frac;
  3866. {
  3867.   unsigned HOST_WIDE_INT ll;
  3868.   unsigned EMUSHORT xi[NI];
  3869.   int j, k;
  3870.  
  3871.   emovi (x, xi);
  3872.   k = (int) xi[E] - (EXONE - 1);
  3873.   if (k <= 0)
  3874.     {
  3875.       /* if exponent <= 0, integer = 0 and argument is fraction */
  3876.       *i = 0L;
  3877.       emovo (xi, frac);
  3878.       return;
  3879.     }
  3880.   if (k > HOST_BITS_PER_WIDE_INT)
  3881.     {
  3882.       /* Long integer overflow: output large integer
  3883.      and correct fraction.
  3884.      Note, the BSD microvax compiler says that ~(0UL)
  3885.      is a syntax error.  */
  3886.       *i = ~(0L);
  3887.       eshift (xi, k);
  3888.       if (extra_warnings)
  3889.     warning ("overflow on truncation to unsigned integer");
  3890.     }
  3891.   else if (k > 16)
  3892.     {
  3893.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  3894.      then shift up by 16's.  */
  3895.       j = k - ((k >> 4) << 4);
  3896.       eshift (xi, j);
  3897.       ll = xi[M];
  3898.       k -= j;
  3899.       do
  3900.     {
  3901.       eshup6 (xi);
  3902.       ll = (ll << 16) | xi[M];
  3903.     }
  3904.       while ((k -= 16) > 0);
  3905.       *i = ll;
  3906.     }
  3907.   else
  3908.     {
  3909.       /* shift not more than 16 bits */
  3910.       eshift (xi, k);
  3911.       *i = (HOST_WIDE_INT) xi[M] & 0xffff;
  3912.     }
  3913.  
  3914.   if (xi[0])  /* A negative value yields unsigned integer 0. */
  3915.     *i = 0L;
  3916.  
  3917.   xi[0] = 0;
  3918.   xi[E] = EXONE - 1;
  3919.   xi[M] = 0;
  3920.   if ((k = enormlz (xi)) > NBITS)
  3921.     ecleaz (xi);
  3922.   else
  3923.     xi[E] -= (unsigned EMUSHORT) k;
  3924.  
  3925.   emovo (xi, frac);
  3926. }
  3927.  
  3928.  
  3929.  
  3930. /* Shift significand area up or down by the number of bits given by SC.  */
  3931.  
  3932. static int 
  3933. eshift (x, sc)
  3934.      unsigned EMUSHORT *x;
  3935.      int sc;
  3936. {
  3937.   unsigned EMUSHORT lost;
  3938.   unsigned EMUSHORT *p;
  3939.  
  3940.   if (sc == 0)
  3941.     return (0);
  3942.  
  3943.   lost = 0;
  3944.   p = x + NI - 1;
  3945.  
  3946.   if (sc < 0)
  3947.     {
  3948.       sc = -sc;
  3949.       while (sc >= 16)
  3950.     {
  3951.       lost |= *p;        /* remember lost bits */
  3952.       eshdn6 (x);
  3953.       sc -= 16;
  3954.     }
  3955.  
  3956.       while (sc >= 8)
  3957.     {
  3958.       lost |= *p & 0xff;
  3959.       eshdn8 (x);
  3960.       sc -= 8;
  3961.     }
  3962.  
  3963.       while (sc > 0)
  3964.     {
  3965.       lost |= *p & 1;
  3966.       eshdn1 (x);
  3967.       sc -= 1;
  3968.     }
  3969.     }
  3970.   else
  3971.     {
  3972.       while (sc >= 16)
  3973.     {
  3974.       eshup6 (x);
  3975.       sc -= 16;
  3976.     }
  3977.  
  3978.       while (sc >= 8)
  3979.     {
  3980.       eshup8 (x);
  3981.       sc -= 8;
  3982.     }
  3983.  
  3984.       while (sc > 0)
  3985.     {
  3986.       eshup1 (x);
  3987.       sc -= 1;
  3988.     }
  3989.     }
  3990.   if (lost)
  3991.     lost = 1;
  3992.   return ((int) lost);
  3993. }
  3994.  
  3995.  
  3996.  
  3997. /* Shift normalize the significand area pointed to by argument.
  3998.    Shift count (up = positive) is returned.  */
  3999.  
  4000. static int 
  4001. enormlz (x)
  4002.      unsigned EMUSHORT x[];
  4003. {
  4004.   register unsigned EMUSHORT *p;
  4005.   int sc;
  4006.  
  4007.   sc = 0;
  4008.   p = &x[M];
  4009.   if (*p != 0)
  4010.     goto normdn;
  4011.   ++p;
  4012.   if (*p & 0x8000)
  4013.     return (0);            /* already normalized */
  4014.   while (*p == 0)
  4015.     {
  4016.       eshup6 (x);
  4017.       sc += 16;
  4018.  
  4019.       /* With guard word, there are NBITS+16 bits available.
  4020.        Return true if all are zero.  */
  4021.       if (sc > NBITS)
  4022.     return (sc);
  4023.     }
  4024.   /* see if high byte is zero */
  4025.   while ((*p & 0xff00) == 0)
  4026.     {
  4027.       eshup8 (x);
  4028.       sc += 8;
  4029.     }
  4030.   /* now shift 1 bit at a time */
  4031.   while ((*p & 0x8000) == 0)
  4032.     {
  4033.       eshup1 (x);
  4034.       sc += 1;
  4035.       if (sc > NBITS)
  4036.     {
  4037.       mtherr ("enormlz", UNDERFLOW);
  4038.       return (sc);
  4039.     }
  4040.     }
  4041.   return (sc);
  4042.  
  4043.   /* Normalize by shifting down out of the high guard word
  4044.      of the significand */
  4045.  normdn:
  4046.  
  4047.   if (*p & 0xff00)
  4048.     {
  4049.       eshdn8 (x);
  4050.       sc -= 8;
  4051.     }
  4052.   while (*p != 0)
  4053.     {
  4054.       eshdn1 (x);
  4055.       sc -= 1;
  4056.  
  4057.       if (sc < -NBITS)
  4058.     {
  4059.       mtherr ("enormlz", OVERFLOW);
  4060.       return (sc);
  4061.     }
  4062.     }
  4063.   return (sc);
  4064. }
  4065.  
  4066.  
  4067.  
  4068.  
  4069. /* Convert e type number to decimal format ASCII string.
  4070.    The constants are for 64 bit precision.  */
  4071.  
  4072. #define NTEN 12
  4073. #define MAXP 4096
  4074.  
  4075. #if LONG_DOUBLE_TYPE_SIZE == 128
  4076. static unsigned EMUSHORT etens[NTEN + 1][NE] =
  4077. {
  4078.   {0x6576, 0x4a92, 0x804a, 0x153f,
  4079.    0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
  4080.   {0x6a32, 0xce52, 0x329a, 0x28ce,
  4081.    0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
  4082.   {0x526c, 0x50ce, 0xf18b, 0x3d28,
  4083.    0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
  4084.   {0x9c66, 0x58f8, 0xbc50, 0x5c54,
  4085.    0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
  4086.   {0x851e, 0xeab7, 0x98fe, 0x901b,
  4087.    0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
  4088.   {0x0235, 0x0137, 0x36b1, 0x336c,
  4089.    0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
  4090.   {0x50f8, 0x25fb, 0xc76b, 0x6b71,
  4091.    0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
  4092.   {0x0000, 0x0000, 0x0000, 0x0000,
  4093.    0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
  4094.   {0x0000, 0x0000, 0x0000, 0x0000,
  4095.    0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
  4096.   {0x0000, 0x0000, 0x0000, 0x0000,
  4097.    0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
  4098.   {0x0000, 0x0000, 0x0000, 0x0000,
  4099.    0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
  4100.   {0x0000, 0x0000, 0x0000, 0x0000,
  4101.    0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
  4102.   {0x0000, 0x0000, 0x0000, 0x0000,
  4103.    0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
  4104. };
  4105.  
  4106. static unsigned EMUSHORT emtens[NTEN + 1][NE] =
  4107. {
  4108.   {0x2030, 0xcffc, 0xa1c3, 0x8123,
  4109.    0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
  4110.   {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
  4111.    0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
  4112.   {0xf53f, 0xf698, 0x6bd3, 0x0158,
  4113.    0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
  4114.   {0xe731, 0x04d4, 0xe3f2, 0xd332,
  4115.    0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
  4116.   {0xa23e, 0x5308, 0xfefb, 0x1155,
  4117.    0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
  4118.   {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
  4119.    0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
  4120.   {0x2a20, 0x6224, 0x47b3, 0x98d7,
  4121.    0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
  4122.   {0x0b5b, 0x4af2, 0xa581, 0x18ed,
  4123.    0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
  4124.   {0xbf71, 0xa9b3, 0x7989, 0xbe68,
  4125.    0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
  4126.   {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
  4127.    0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
  4128.   {0xc155, 0xa4a8, 0x404e, 0x6113,
  4129.    0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
  4130.   {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
  4131.    0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
  4132.   {0xcccd, 0xcccc, 0xcccc, 0xcccc,
  4133.    0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
  4134. };
  4135. #else
  4136. /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
  4137. static unsigned EMUSHORT etens[NTEN + 1][NE] =
  4138. {
  4139.   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
  4140.   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
  4141.   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
  4142.   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
  4143.   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
  4144.   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
  4145.   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
  4146.   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
  4147.   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
  4148.   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
  4149.   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
  4150.   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
  4151.   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
  4152. };
  4153.  
  4154. static unsigned EMUSHORT emtens[NTEN + 1][NE] =
  4155. {
  4156.   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
  4157.   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
  4158.   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
  4159.   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
  4160.   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
  4161.   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
  4162.   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
  4163.   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
  4164.   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
  4165.   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
  4166.   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
  4167.   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
  4168.   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
  4169. };
  4170. #endif
  4171.  
  4172. static void 
  4173. e24toasc (x, string, ndigs)
  4174.      unsigned EMUSHORT x[];
  4175.      char *string;
  4176.      int ndigs;
  4177. {
  4178.   unsigned EMUSHORT w[NI];
  4179.  
  4180.   e24toe (x, w);
  4181.   etoasc (w, string, ndigs);
  4182. }
  4183.  
  4184.  
  4185. static void 
  4186. e53toasc (x, string, ndigs)
  4187.      unsigned EMUSHORT x[];
  4188.      char *string;
  4189.      int ndigs;
  4190. {
  4191.   unsigned EMUSHORT w[NI];
  4192.  
  4193.   e53toe (x, w);
  4194.   etoasc (w, string, ndigs);
  4195. }
  4196.  
  4197.  
  4198. static void 
  4199. e64toasc (x, string, ndigs)
  4200.      unsigned EMUSHORT x[];
  4201.      char *string;
  4202.      int ndigs;
  4203. {
  4204.   unsigned EMUSHORT w[NI];
  4205.  
  4206.   e64toe (x, w);
  4207.   etoasc (w, string, ndigs);
  4208. }
  4209.  
  4210. static void 
  4211. e113toasc (x, string, ndigs)
  4212.      unsigned EMUSHORT x[];
  4213.      char *string;
  4214.      int ndigs;
  4215. {
  4216.   unsigned EMUSHORT w[NI];
  4217.  
  4218.   e113toe (x, w);
  4219.   etoasc (w, string, ndigs);
  4220. }
  4221.  
  4222.  
  4223. static char wstring[80];    /* working storage for ASCII output */
  4224.  
  4225. static void 
  4226. etoasc (x, string, ndigs)
  4227.      unsigned EMUSHORT x[];
  4228.      char *string;
  4229.      int ndigs;
  4230. {
  4231.   EMUSHORT digit;
  4232.   unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
  4233.   unsigned EMUSHORT *p, *r, *ten;
  4234.   unsigned EMUSHORT sign;
  4235.   int i, j, k, expon, rndsav;
  4236.   char *s, *ss;
  4237.   unsigned EMUSHORT m;
  4238.  
  4239.  
  4240.   rndsav = rndprc;
  4241.   ss = string;
  4242.   s = wstring;
  4243.   *ss = '\0';
  4244.   *s = '\0';
  4245. #ifdef NANS
  4246.   if (eisnan (x))
  4247.     {
  4248.       sprintf (wstring, " NaN ");
  4249.       goto bxit;
  4250.     }
  4251. #endif
  4252.   rndprc = NBITS;        /* set to full precision */
  4253.   emov (x, y);            /* retain external format */
  4254.   if (y[NE - 1] & 0x8000)
  4255.     {
  4256.       sign = 0xffff;
  4257.       y[NE - 1] &= 0x7fff;
  4258.     }
  4259.   else
  4260.     {
  4261.       sign = 0;
  4262.     }
  4263.   expon = 0;
  4264.   ten = &etens[NTEN][0];
  4265.   emov (eone, t);
  4266.   /* Test for zero exponent */
  4267.   if (y[NE - 1] == 0)
  4268.     {
  4269.       for (k = 0; k < NE - 1; k++)
  4270.     {
  4271.       if (y[k] != 0)
  4272.         goto tnzro;        /* denormalized number */
  4273.     }
  4274.       goto isone;        /* legal all zeros */
  4275.     }
  4276.  tnzro:
  4277.  
  4278.   /* Test for infinity. */
  4279.   if (y[NE - 1] == 0x7fff)
  4280.     {
  4281.       if (sign)
  4282.     sprintf (wstring, " -Infinity ");
  4283.       else
  4284.     sprintf (wstring, " Infinity ");
  4285.       goto bxit;
  4286.     }
  4287.  
  4288.   /* Test for exponent nonzero but significand denormalized.
  4289.    * This is an error condition.
  4290.    */
  4291.   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
  4292.     {
  4293.       mtherr ("etoasc", DOMAIN);
  4294.       sprintf (wstring, "NaN");
  4295.       goto bxit;
  4296.     }
  4297.  
  4298.   /* Compare to 1.0 */
  4299.   i = ecmp (eone, y);
  4300.   if (i == 0)
  4301.     goto isone;
  4302.  
  4303.   if (i == -2)
  4304.     abort ();
  4305.  
  4306.   if (i < 0)
  4307.     {                /* Number is greater than 1 */
  4308.       /* Convert significand to an integer and strip trailing decimal zeros. */
  4309.       emov (y, u);
  4310.       u[NE - 1] = EXONE + NBITS - 1;
  4311.  
  4312.       p = &etens[NTEN - 4][0];
  4313.       m = 16;
  4314.       do
  4315.     {
  4316.       ediv (p, u, t);
  4317.       efloor (t, w);
  4318.       for (j = 0; j < NE - 1; j++)
  4319.         {
  4320.           if (t[j] != w[j])
  4321.         goto noint;
  4322.         }
  4323.       emov (t, u);
  4324.       expon += (int) m;
  4325.     noint:
  4326.       p += NE;
  4327.       m >>= 1;
  4328.     }
  4329.       while (m != 0);
  4330.  
  4331.       /* Rescale from integer significand */
  4332.       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
  4333.       emov (u, y);
  4334.       /* Find power of 10 */
  4335.       emov (eone, t);
  4336.       m = MAXP;
  4337.       p = &etens[0][0];
  4338.       /* An unordered compare result shouldn't happen here. */
  4339.       while (ecmp (ten, u) <= 0)
  4340.     {
  4341.       if (ecmp (p, u) <= 0)
  4342.         {
  4343.           ediv (p, u, u);
  4344.           emul (p, t, t);
  4345.           expon += (int) m;
  4346.         }
  4347.       m >>= 1;
  4348.       if (m == 0)
  4349.         break;
  4350.       p += NE;
  4351.     }
  4352.     }
  4353.   else
  4354.     {                /* Number is less than 1.0 */
  4355.       /* Pad significand with trailing decimal zeros. */
  4356.       if (y[NE - 1] == 0)
  4357.     {
  4358.       while ((y[NE - 2] & 0x8000) == 0)
  4359.         {
  4360.           emul (ten, y, y);
  4361.           expon -= 1;
  4362.         }
  4363.     }
  4364.       else
  4365.     {
  4366.       emovi (y, w);
  4367.       for (i = 0; i < NDEC + 1; i++)
  4368.         {
  4369.           if ((w[NI - 1] & 0x7) != 0)
  4370.         break;
  4371.           /* multiply by 10 */
  4372.           emovz (w, u);
  4373.           eshdn1 (u);
  4374.           eshdn1 (u);
  4375.           eaddm (w, u);
  4376.           u[1] += 3;
  4377.           while (u[2] != 0)
  4378.         {
  4379.           eshdn1 (u);
  4380.           u[1] += 1;
  4381.         }
  4382.           if (u[NI - 1] != 0)
  4383.         break;
  4384.           if (eone[NE - 1] <= u[1])
  4385.         break;
  4386.           emovz (u, w);
  4387.           expon -= 1;
  4388.         }
  4389.       emovo (w, y);
  4390.     }
  4391.       k = -MAXP;
  4392.       p = &emtens[0][0];
  4393.       r = &etens[0][0];
  4394.       emov (y, w);
  4395.       emov (eone, t);
  4396.       while (ecmp (eone, w) > 0)
  4397.     {
  4398.       if (ecmp (p, w) >= 0)
  4399.         {
  4400.           emul (r, w, w);
  4401.           emul (r, t, t);
  4402.           expon += k;
  4403.         }
  4404.       k /= 2;
  4405.       if (k == 0)
  4406.         break;
  4407.       p += NE;
  4408.       r += NE;
  4409.     }
  4410.       ediv (t, eone, t);
  4411.     }
  4412.  isone:
  4413.   /* Find the first (leading) digit. */
  4414.   emovi (t, w);
  4415.   emovz (w, t);
  4416.   emovi (y, w);
  4417.   emovz (w, y);
  4418.   eiremain (t, y);
  4419.   digit = equot[NI - 1];
  4420.   while ((digit == 0) && (ecmp (y, ezero) != 0))
  4421.     {
  4422.       eshup1 (y);
  4423.       emovz (y, u);
  4424.       eshup1 (u);
  4425.       eshup1 (u);
  4426.       eaddm (u, y);
  4427.       eiremain (t, y);
  4428.       digit = equot[NI - 1];
  4429.       expon -= 1;
  4430.     }
  4431.   s = wstring;
  4432.   if (sign)
  4433.     *s++ = '-';
  4434.   else
  4435.     *s++ = ' ';
  4436.   /* Examine number of digits requested by caller. */
  4437.   if (ndigs < 0)
  4438.     ndigs = 0;
  4439.   if (ndigs > NDEC)
  4440.     ndigs = NDEC;
  4441.   if (digit == 10)
  4442.     {
  4443.       *s++ = '1';
  4444.       *s++ = '.';
  4445.       if (ndigs > 0)
  4446.     {
  4447.       *s++ = '0';
  4448.       ndigs -= 1;
  4449.     }
  4450.       expon += 1;
  4451.     }
  4452.   else
  4453.     {
  4454.       *s++ = (char)digit + '0';
  4455.       *s++ = '.';
  4456.     }
  4457.   /* Generate digits after the decimal point. */
  4458.   for (k = 0; k <= ndigs; k++)
  4459.     {
  4460.       /* multiply current number by 10, without normalizing */
  4461.       eshup1 (y);
  4462.       emovz (y, u);
  4463.       eshup1 (u);
  4464.       eshup1 (u);
  4465.       eaddm (u, y);
  4466.       eiremain (t, y);
  4467.       *s++ = (char) equot[NI - 1] + '0';
  4468.     }
  4469.   digit = equot[NI - 1];
  4470.   --s;
  4471.   ss = s;
  4472.   /* round off the ASCII string */
  4473.   if (digit > 4)
  4474.     {
  4475.       /* Test for critical rounding case in ASCII output. */
  4476.       if (digit == 5)
  4477.     {
  4478.       emovo (y, t);
  4479.       if (ecmp (t, ezero) != 0)
  4480.         goto roun;        /* round to nearest */
  4481.       if ((*(s - 1) & 1) == 0)
  4482.         goto doexp;        /* round to even */
  4483.     }
  4484.       /* Round up and propagate carry-outs */
  4485.     roun:
  4486.       --s;
  4487.       k = *s & 0x7f;
  4488.       /* Carry out to most significant digit? */
  4489.       if (k == '.')
  4490.     {
  4491.       --s;
  4492.       k = *s;
  4493.       k += 1;
  4494.       *s = (char) k;
  4495.       /* Most significant digit carries to 10? */
  4496.       if (k > '9')
  4497.         {
  4498.           expon += 1;
  4499.           *s = '1';
  4500.         }
  4501.       goto doexp;
  4502.     }
  4503.       /* Round up and carry out from less significant digits */
  4504.       k += 1;
  4505.       *s = (char) k;
  4506.       if (k > '9')
  4507.     {
  4508.       *s = '0';
  4509.       goto roun;
  4510.     }
  4511.     }
  4512.  doexp:
  4513.   /*
  4514.      if (expon >= 0)
  4515.      sprintf (ss, "e+%d", expon);
  4516.      else
  4517.      sprintf (ss, "e%d", expon);
  4518.      */
  4519.   sprintf (ss, "e%d", expon);
  4520.  bxit:
  4521.   rndprc = rndsav;
  4522.   /* copy out the working string */
  4523.   s = string;
  4524.   ss = wstring;
  4525.   while (*ss == ' ')        /* strip possible leading space */
  4526.     ++ss;
  4527.   while ((*s++ = *ss++) != '\0')
  4528.     ;
  4529. }
  4530.  
  4531.  
  4532. /* Convert ASCII string to quadruple precision floating point
  4533.  
  4534.    Numeric input is free field decimal number with max of 15 digits with or
  4535.    without decimal point entered as ASCII from teletype.  Entering E after
  4536.    the number followed by a second number causes the second number to be
  4537.    interpreted as a power of 10 to be multiplied by the first number
  4538.    (i.e., "scientific" notation).  */
  4539.  
  4540. /* ASCII to single */
  4541.  
  4542. static void 
  4543. asctoe24 (s, y)
  4544.      char *s;
  4545.      unsigned EMUSHORT *y;
  4546. {
  4547.   asctoeg (s, y, 24);
  4548. }
  4549.  
  4550.  
  4551. /* ASCII to double */
  4552.  
  4553. static void 
  4554. asctoe53 (s, y)
  4555.      char *s;
  4556.      unsigned EMUSHORT *y;
  4557. {
  4558. #if defined(DEC) || defined(IBM)
  4559.   asctoeg (s, y, 56);
  4560. #else
  4561.   asctoeg (s, y, 53);
  4562. #endif
  4563. }
  4564.  
  4565.  
  4566. /* ASCII to long double */
  4567.  
  4568. static void 
  4569. asctoe64 (s, y)
  4570.      char *s;
  4571.      unsigned EMUSHORT *y;
  4572. {
  4573.   asctoeg (s, y, 64);
  4574. }
  4575.  
  4576. /* ASCII to 128-bit long double */
  4577.  
  4578. static void 
  4579. asctoe113 (s, y)
  4580.      char *s;
  4581.      unsigned EMUSHORT *y;
  4582. {
  4583.   asctoeg (s, y, 113);
  4584. }
  4585.  
  4586. /* ASCII to super double */
  4587.  
  4588. static void 
  4589. asctoe (s, y)
  4590.      char *s;
  4591.      unsigned EMUSHORT *y;
  4592. {
  4593.   asctoeg (s, y, NBITS);
  4594. }
  4595.  
  4596.  
  4597. /* ASCII to e type, with specified rounding precision = oprec. */
  4598.  
  4599. static void 
  4600. asctoeg (ss, y, oprec)
  4601.      char *ss;
  4602.      unsigned EMUSHORT *y;
  4603.      int oprec;
  4604. {
  4605.   unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
  4606.   int esign, decflg, sgnflg, nexp, exp, prec, lost;
  4607.   int k, trail, c, rndsav;
  4608.   EMULONG lexp;
  4609.   unsigned EMUSHORT nsign, *p;
  4610.   char *sp, *s, *lstr;
  4611.  
  4612.   /* Copy the input string. */
  4613.   lstr = (char *) alloca (strlen (ss) + 1);
  4614.   s = ss;
  4615.   while (*s == ' ')        /* skip leading spaces */
  4616.     ++s;
  4617.   sp = lstr;
  4618.   while ((*sp++ = *s++) != '\0')
  4619.     ;
  4620.   s = lstr;
  4621.  
  4622.   rndsav = rndprc;
  4623.   rndprc = NBITS;        /* Set to full precision */
  4624.   lost = 0;
  4625.   nsign = 0;
  4626.   decflg = 0;
  4627.   sgnflg = 0;
  4628.   nexp = 0;
  4629.   exp = 0;
  4630.   prec = 0;
  4631.   ecleaz (yy);
  4632.   trail = 0;
  4633.  
  4634.  nxtcom:
  4635.   k = *s - '0';
  4636.   if ((k >= 0) && (k <= 9))
  4637.     {
  4638.       /* Ignore leading zeros */
  4639.       if ((prec == 0) && (decflg == 0) && (k == 0))
  4640.     goto donchr;
  4641.       /* Identify and strip trailing zeros after the decimal point. */
  4642.       if ((trail == 0) && (decflg != 0))
  4643.     {
  4644.       sp = s;
  4645.       while ((*sp >= '0') && (*sp <= '9'))
  4646.         ++sp;
  4647.       /* Check for syntax error */
  4648.       c = *sp & 0x7f;
  4649.       if ((c != 'e') && (c != 'E') && (c != '\0')
  4650.           && (c != '\n') && (c != '\r') && (c != ' ')
  4651.           && (c != ','))
  4652.         goto error;
  4653.       --sp;
  4654.       while (*sp == '0')
  4655.         *sp-- = 'z';
  4656.       trail = 1;
  4657.       if (*s == 'z')
  4658.         goto donchr;
  4659.     }
  4660.  
  4661.       /* If enough digits were given to more than fill up the yy register,
  4662.      continuing until overflow into the high guard word yy[2]
  4663.      guarantees that there will be a roundoff bit at the top
  4664.      of the low guard word after normalization.  */
  4665.  
  4666.       if (yy[2] == 0)
  4667.     {
  4668.       if (decflg)
  4669.         nexp += 1;        /* count digits after decimal point */
  4670.       eshup1 (yy);        /* multiply current number by 10 */
  4671.       emovz (yy, xt);
  4672.       eshup1 (xt);
  4673.       eshup1 (xt);
  4674.       eaddm (xt, yy);
  4675.       ecleaz (xt);
  4676.       xt[NI - 2] = (unsigned EMUSHORT) k;
  4677.       eaddm (xt, yy);
  4678.     }
  4679.       else
  4680.     {
  4681.       /* Mark any lost non-zero digit.  */
  4682.       lost |= k;
  4683.       /* Count lost digits before the decimal point.  */
  4684.       if (decflg == 0)
  4685.         nexp -= 1;
  4686.     }
  4687.       prec += 1;
  4688.       goto donchr;
  4689.     }
  4690.  
  4691.   switch (*s)
  4692.     {
  4693.     case 'z':
  4694.       break;
  4695.     case 'E':
  4696.     case 'e':
  4697.       goto expnt;
  4698.     case '.':            /* decimal point */
  4699.       if (decflg)
  4700.     goto error;
  4701.       ++decflg;
  4702.       break;
  4703.     case '-':
  4704.       nsign = 0xffff;
  4705.       if (sgnflg)
  4706.     goto error;
  4707.       ++sgnflg;
  4708.       break;
  4709.     case '+':
  4710.       if (sgnflg)
  4711.     goto error;
  4712.       ++sgnflg;
  4713.       break;
  4714.     case ',':
  4715.     case ' ':
  4716.     case '\0':
  4717.     case '\n':
  4718.     case '\r':
  4719.       goto daldone;
  4720.     case 'i':
  4721.     case 'I':
  4722.       goto infinite;
  4723.     default:
  4724.     error:
  4725. #ifdef NANS
  4726.       einan (yy);
  4727. #else
  4728.       mtherr ("asctoe", DOMAIN);
  4729.       eclear (yy);
  4730. #endif
  4731.       goto aexit;
  4732.     }
  4733.  donchr:
  4734.   ++s;
  4735.   goto nxtcom;
  4736.  
  4737.   /* Exponent interpretation */
  4738.  expnt:
  4739.  
  4740.   esign = 1;
  4741.   exp = 0;
  4742.   ++s;
  4743.   /* check for + or - */
  4744.   if (*s == '-')
  4745.     {
  4746.       esign = -1;
  4747.       ++s;
  4748.     }
  4749.   if (*s == '+')
  4750.     ++s;
  4751.   while ((*s >= '0') && (*s <= '9'))
  4752.     {
  4753.       exp *= 10;
  4754.       exp += *s++ - '0';
  4755.       if (exp > -(MINDECEXP))
  4756.     {
  4757.       if (esign < 0)
  4758.         goto zero;
  4759.       else
  4760.         goto infinite;
  4761.     }
  4762.     }
  4763.   if (esign < 0)
  4764.     exp = -exp;
  4765.   if (exp > MAXDECEXP)
  4766.     {
  4767.  infinite:
  4768.       ecleaz (yy);
  4769.       yy[E] = 0x7fff;        /* infinity */
  4770.       goto aexit;
  4771.     }
  4772.   if (exp < MINDECEXP)
  4773.     {
  4774.  zero:
  4775.       ecleaz (yy);
  4776.       goto aexit;
  4777.     }
  4778.  
  4779.  daldone:
  4780.   nexp = exp - nexp;
  4781.   /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
  4782.   while ((nexp > 0) && (yy[2] == 0))
  4783.     {
  4784.       emovz (yy, xt);
  4785.       eshup1 (xt);
  4786.       eshup1 (xt);
  4787.       eaddm (yy, xt);
  4788.       eshup1 (xt);
  4789.       if (xt[2] != 0)
  4790.     break;
  4791.       nexp -= 1;
  4792.       emovz (xt, yy);
  4793.     }
  4794.   if ((k = enormlz (yy)) > NBITS)
  4795.     {
  4796.       ecleaz (yy);
  4797.       goto aexit;
  4798.     }
  4799.   lexp = (EXONE - 1 + NBITS) - k;
  4800.   emdnorm (yy, lost, 0, lexp, 64);
  4801.  
  4802.   /* Convert to external format:
  4803.  
  4804.      Multiply by 10**nexp.  If precision is 64 bits,
  4805.      the maximum relative error incurred in forming 10**n
  4806.      for 0 <= n <= 324 is 8.2e-20, at 10**180.
  4807.      For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
  4808.      For 0 >= n >= -999, it is -1.55e-19 at 10**-435.  */
  4809.  
  4810.   lexp = yy[E];
  4811.   if (nexp == 0)
  4812.     {
  4813.       k = 0;
  4814.       goto expdon;
  4815.     }
  4816.   esign = 1;
  4817.   if (nexp < 0)
  4818.     {
  4819.       nexp = -nexp;
  4820.       esign = -1;
  4821.       if (nexp > 4096)
  4822.     {
  4823.       /* Punt.  Can't handle this without 2 divides. */
  4824.       emovi (etens[0], tt);
  4825.       lexp -= tt[E];
  4826.       k = edivm (tt, yy);
  4827.       lexp += EXONE;
  4828.       nexp -= 4096;
  4829.     }
  4830.     }
  4831.   p = &etens[NTEN][0];
  4832.   emov (eone, xt);
  4833.   exp = 1;
  4834.   do
  4835.     {
  4836.       if (exp & nexp)
  4837.     emul (p, xt, xt);
  4838.       p -= NE;
  4839.       exp = exp + exp;
  4840.     }
  4841.   while (exp <= MAXP);
  4842.  
  4843.   emovi (xt, tt);
  4844.   if (esign < 0)
  4845.     {
  4846.       lexp -= tt[E];
  4847.       k = edivm (tt, yy);
  4848.       lexp += EXONE;
  4849.     }
  4850.   else
  4851.     {
  4852.       lexp += tt[E];
  4853.       k = emulm (tt, yy);
  4854.       lexp -= EXONE - 1;
  4855.     }
  4856.  
  4857.  expdon:
  4858.  
  4859.   /* Round and convert directly to the destination type */
  4860.   if (oprec == 53)
  4861.     lexp -= EXONE - 0x3ff;
  4862. #ifdef IBM
  4863.   else if (oprec == 24 || oprec == 56)
  4864.     lexp -= EXONE - (0x41 << 2);
  4865. #else
  4866.   else if (oprec == 24)
  4867.     lexp -= EXONE - 0177;
  4868. #endif
  4869. #ifdef DEC
  4870.   else if (oprec == 56)
  4871.     lexp -= EXONE - 0201;
  4872. #endif
  4873.   rndprc = oprec;
  4874.   emdnorm (yy, k, 0, lexp, 64);
  4875.  
  4876.  aexit:
  4877.  
  4878.   rndprc = rndsav;
  4879.   yy[0] = nsign;
  4880.   switch (oprec)
  4881.     {
  4882. #ifdef DEC
  4883.     case 56:
  4884.       todec (yy, y);        /* see etodec.c */
  4885.       break;
  4886. #endif
  4887. #ifdef IBM
  4888.     case 56:
  4889.       toibm (yy, y, DFmode);
  4890.       break;
  4891. #endif
  4892.     case 53:
  4893.       toe53 (yy, y);
  4894.       break;
  4895.     case 24:
  4896.       toe24 (yy, y);
  4897.       break;
  4898.     case 64:
  4899.       toe64 (yy, y);
  4900.       break;
  4901.     case 113:
  4902.       toe113 (yy, y);
  4903.       break;
  4904.     case NBITS:
  4905.       emovo (yy, y);
  4906.       break;
  4907.     }
  4908. }
  4909.  
  4910.  
  4911.  
  4912. /* y = largest integer not greater than x (truncated toward minus infinity)  */
  4913.  
  4914. static unsigned EMUSHORT bmask[] =
  4915. {
  4916.   0xffff,
  4917.   0xfffe,
  4918.   0xfffc,
  4919.   0xfff8,
  4920.   0xfff0,
  4921.   0xffe0,
  4922.   0xffc0,
  4923.   0xff80,
  4924.   0xff00,
  4925.   0xfe00,
  4926.   0xfc00,
  4927.   0xf800,
  4928.   0xf000,
  4929.   0xe000,
  4930.   0xc000,
  4931.   0x8000,
  4932.   0x0000,
  4933. };
  4934.  
  4935. static void 
  4936. efloor (x, y)
  4937.      unsigned EMUSHORT x[], y[];
  4938. {
  4939.   register unsigned EMUSHORT *p;
  4940.   int e, expon, i;
  4941.   unsigned EMUSHORT f[NE];
  4942.  
  4943.   emov (x, f);            /* leave in external format */
  4944.   expon = (int) f[NE - 1];
  4945.   e = (expon & 0x7fff) - (EXONE - 1);
  4946.   if (e <= 0)
  4947.     {
  4948.       eclear (y);
  4949.       goto isitneg;
  4950.     }
  4951.   /* number of bits to clear out */
  4952.   e = NBITS - e;
  4953.   emov (f, y);
  4954.   if (e <= 0)
  4955.     return;
  4956.  
  4957.   p = &y[0];
  4958.   while (e >= 16)
  4959.     {
  4960.       *p++ = 0;
  4961.       e -= 16;
  4962.     }
  4963.   /* clear the remaining bits */
  4964.   *p &= bmask[e];
  4965.   /* truncate negatives toward minus infinity */
  4966.  isitneg:
  4967.  
  4968.   if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
  4969.     {
  4970.       for (i = 0; i < NE - 1; i++)
  4971.     {
  4972.       if (f[i] != y[i])
  4973.         {
  4974.           esub (eone, y, y);
  4975.           break;
  4976.         }
  4977.     }
  4978.     }
  4979. }
  4980.  
  4981.  
  4982. /* Returns s and exp such that  s * 2**exp = x and .5 <= s < 1.
  4983.    For example, 1.1 = 0.55 * 2**1
  4984.    Handles denormalized numbers properly using long integer exp.  */
  4985.  
  4986. static void 
  4987. efrexp (x, exp, s)
  4988.      unsigned EMUSHORT x[];
  4989.      int *exp;
  4990.      unsigned EMUSHORT s[];
  4991. {
  4992.   unsigned EMUSHORT xi[NI];
  4993.   EMULONG li;
  4994.  
  4995.   emovi (x, xi);
  4996.   li = (EMULONG) ((EMUSHORT) xi[1]);
  4997.  
  4998.   if (li == 0)
  4999.     {
  5000.       li -= enormlz (xi);
  5001.     }
  5002.   xi[1] = 0x3ffe;
  5003.   emovo (xi, s);
  5004.   *exp = (int) (li - 0x3ffe);
  5005. }
  5006.  
  5007.  
  5008.  
  5009. /* Return y = x * 2**pwr2.  */
  5010.  
  5011. static void 
  5012. eldexp (x, pwr2, y)
  5013.      unsigned EMUSHORT x[];
  5014.      int pwr2;
  5015.      unsigned EMUSHORT y[];
  5016. {
  5017.   unsigned EMUSHORT xi[NI];
  5018.   EMULONG li;
  5019.   int i;
  5020.  
  5021.   emovi (x, xi);
  5022.   li = xi[1];
  5023.   li += pwr2;
  5024.   i = 0;
  5025.   emdnorm (xi, i, i, li, 64);
  5026.   emovo (xi, y);
  5027. }
  5028.  
  5029.  
  5030. /* c = remainder after dividing b by a
  5031.    Least significant integer quotient bits left in equot[].  */
  5032.  
  5033. static void 
  5034. eremain (a, b, c)
  5035.      unsigned EMUSHORT a[], b[], c[];
  5036. {
  5037.   unsigned EMUSHORT den[NI], num[NI];
  5038.  
  5039. #ifdef NANS
  5040.   if (eisinf (b)
  5041.       || (ecmp (a, ezero) == 0)
  5042.       || eisnan (a)
  5043.       || eisnan (b))
  5044.     {
  5045.       enan (c, 0);
  5046.       return;
  5047.     }
  5048. #endif
  5049.   if (ecmp (a, ezero) == 0)
  5050.     {
  5051.       mtherr ("eremain", SING);
  5052.       eclear (c);
  5053.       return;
  5054.     }
  5055.   emovi (a, den);
  5056.   emovi (b, num);
  5057.   eiremain (den, num);
  5058.   /* Sign of remainder = sign of quotient */
  5059.   if (a[0] == b[0])
  5060.     num[0] = 0;
  5061.   else
  5062.     num[0] = 0xffff;
  5063.   emovo (num, c);
  5064. }
  5065.  
  5066. static void 
  5067. eiremain (den, num)
  5068.      unsigned EMUSHORT den[], num[];
  5069. {
  5070.   EMULONG ld, ln;
  5071.   unsigned EMUSHORT j;
  5072.  
  5073.   ld = den[E];
  5074.   ld -= enormlz (den);
  5075.   ln = num[E];
  5076.   ln -= enormlz (num);
  5077.   ecleaz (equot);
  5078.   while (ln >= ld)
  5079.     {
  5080.       if (ecmpm (den, num) <= 0)
  5081.     {
  5082.       esubm (den, num);
  5083.       j = 1;
  5084.     }
  5085.       else
  5086.     {
  5087.       j = 0;
  5088.     }
  5089.       eshup1 (equot);
  5090.       equot[NI - 1] |= j;
  5091.       eshup1 (num);
  5092.       ln -= 1;
  5093.     }
  5094.   emdnorm (num, 0, 0, ln, 0);
  5095. }
  5096.  
  5097. /* This routine may be called to report one of the following
  5098.    error conditions (in the include file mconf.h).
  5099.  
  5100.     Mnemonic        Value          Significance
  5101.  
  5102.      DOMAIN            1       argument domain error
  5103.      SING              2       function singularity
  5104.      OVERFLOW          3       overflow range error
  5105.      UNDERFLOW         4       underflow range error
  5106.      TLOSS             5       total loss of precision
  5107.      PLOSS             6       partial loss of precision
  5108.      INVALID           7       NaN - producing operation
  5109.      EDOM             33       Unix domain error code
  5110.      ERANGE           34       Unix range error code
  5111.  
  5112.   The default version of the file prints the function name,
  5113.   passed to it by the pointer fctnam, followed by the
  5114.   error condition.  The display is directed to the standard
  5115.   output device.  The routine then returns to the calling
  5116.   program.  Users may wish to modify the program to abort by
  5117.   calling exit under severe error conditions such as domain
  5118.   errors.
  5119.  
  5120.   Since all error conditions pass control to this function,
  5121.   the display may be easily changed, eliminated, or directed
  5122.   to an error logging device. */
  5123.  
  5124. /* Note: the order of appearance of the following messages is bound to the
  5125.    error codes defined above.  */
  5126.  
  5127. #define NMSGS 8
  5128. static char *ermsg[NMSGS] =
  5129. {
  5130.   "unknown",            /* error code 0 */
  5131.   "domain",            /* error code 1 */
  5132.   "singularity",        /* et seq.      */
  5133.   "overflow",
  5134.   "underflow",
  5135.   "total loss of precision",
  5136.   "partial loss of precision",
  5137.   "invalid operation"
  5138. };
  5139.  
  5140. int merror = 0;
  5141. extern int merror;
  5142.  
  5143. static void 
  5144. mtherr (name, code)
  5145.      char *name;
  5146.      int code;
  5147. {
  5148.   char errstr[80];
  5149.  
  5150.   /* Display string passed by calling program, which is supposed to be the
  5151.      name of the function in which the error occurred.
  5152.  
  5153.      Display error message defined by the code argument.  */
  5154.  
  5155.   if ((code <= 0) || (code >= NMSGS))
  5156.     code = 0;
  5157.   sprintf (errstr, " %s %s error", name, ermsg[code]);
  5158.   if (extra_warnings)
  5159.     warning (errstr);
  5160.   /* Set global error message word */
  5161.   merror = code + 1;
  5162. }
  5163.  
  5164. #ifdef DEC
  5165. /* Convert DEC double precision to e type.  */
  5166.  
  5167. static void 
  5168. dectoe (d, e)
  5169.      unsigned EMUSHORT *d;
  5170.      unsigned EMUSHORT *e;
  5171. {
  5172.   unsigned EMUSHORT y[NI];
  5173.   register unsigned EMUSHORT r, *p;
  5174.  
  5175.   ecleaz (y);            /* start with a zero */
  5176.   p = y;            /* point to our number */
  5177.   r = *d;            /* get DEC exponent word */
  5178.   if (*d & (unsigned int) 0x8000)
  5179.     *p = 0xffff;        /* fill in our sign */
  5180.   ++p;                /* bump pointer to our exponent word */
  5181.   r &= 0x7fff;            /* strip the sign bit */
  5182.   if (r == 0)            /* answer = 0 if high order DEC word = 0 */
  5183.     goto done;
  5184.  
  5185.  
  5186.   r >>= 7;            /* shift exponent word down 7 bits */
  5187.   r += EXONE - 0201;        /* subtract DEC exponent offset */
  5188.   /* add our e type exponent offset */
  5189.   *p++ = r;            /* to form our exponent */
  5190.  
  5191.   r = *d++;            /* now do the high order mantissa */
  5192.   r &= 0177;            /* strip off the DEC exponent and sign bits */
  5193.   r |= 0200;            /* the DEC understood high order mantissa bit */
  5194.   *p++ = r;            /* put result in our high guard word */
  5195.  
  5196.   *p++ = *d++;            /* fill in the rest of our mantissa */
  5197.   *p++ = *d++;
  5198.   *p = *d;
  5199.  
  5200.   eshdn8 (y);            /* shift our mantissa down 8 bits */
  5201.  done:
  5202.   emovo (y, e);
  5203. }
  5204.  
  5205.  
  5206.  
  5207. /*
  5208. ;    convert e type to DEC double precision
  5209. ;    double d;
  5210. ;    EMUSHORT e[NE];
  5211. ;    etodec (e, &d);
  5212. */
  5213.  
  5214. static void 
  5215. etodec (x, d)
  5216.      unsigned EMUSHORT *x, *d;
  5217. {
  5218.   unsigned EMUSHORT xi[NI];
  5219.   EMULONG exp;
  5220.   int rndsav;
  5221.  
  5222.   emovi (x, xi);
  5223.   exp = (EMULONG) xi[E] - (EXONE - 0201);    /* adjust exponent for offsets */
  5224. /* round off to nearest or even */
  5225.   rndsav = rndprc;
  5226.   rndprc = 56;
  5227.   emdnorm (xi, 0, 0, exp, 64);
  5228.   rndprc = rndsav;
  5229.   todec (xi, d);
  5230. }
  5231.  
  5232. static void 
  5233. todec (x, y)
  5234.      unsigned EMUSHORT *x, *y;
  5235. {
  5236.   unsigned EMUSHORT i;
  5237.   unsigned EMUSHORT *p;
  5238.  
  5239.   p = x;
  5240.   *y = 0;
  5241.   if (*p++)
  5242.     *y = 0100000;
  5243.   i = *p++;
  5244.   if (i == 0)
  5245.     {
  5246.       *y++ = 0;
  5247.       *y++ = 0;
  5248.       *y++ = 0;
  5249.       *y++ = 0;
  5250.       return;
  5251.     }
  5252.   if (i > 0377)
  5253.     {
  5254.       *y++ |= 077777;
  5255.       *y++ = 0xffff;
  5256.       *y++ = 0xffff;
  5257.       *y++ = 0xffff;
  5258. #ifdef ERANGE
  5259.       errno = ERANGE;
  5260. #endif
  5261.       return;
  5262.     }
  5263.   i &= 0377;
  5264.   i <<= 7;
  5265.   eshup8 (x);
  5266.   x[M] &= 0177;
  5267.   i |= x[M];
  5268.   *y++ |= i;
  5269.   *y++ = x[M + 1];
  5270.   *y++ = x[M + 2];
  5271.   *y++ = x[M + 3];
  5272. }
  5273. #endif /* DEC */
  5274.  
  5275. #ifdef IBM
  5276. /* Convert IBM single/double precision to e type.  */
  5277.  
  5278. static void 
  5279. ibmtoe (d, e, mode)
  5280.      unsigned EMUSHORT *d;
  5281.      unsigned EMUSHORT *e;
  5282.      enum machine_mode mode;
  5283. {
  5284.   unsigned EMUSHORT y[NI];
  5285.   register unsigned EMUSHORT r, *p;
  5286.   int rndsav;
  5287.  
  5288.   ecleaz (y);            /* start with a zero */
  5289.   p = y;            /* point to our number */
  5290.   r = *d;            /* get IBM exponent word */
  5291.   if (*d & (unsigned int) 0x8000)
  5292.     *p = 0xffff;        /* fill in our sign */
  5293.   ++p;                /* bump pointer to our exponent word */
  5294.   r &= 0x7f00;            /* strip the sign bit */
  5295.   r >>= 6;            /* shift exponent word down 6 bits */
  5296.                 /* in fact shift by 8 right and 2 left */
  5297.   r += EXONE - (0x41 << 2);    /* subtract IBM exponent offset */
  5298.                   /* add our e type exponent offset */
  5299.   *p++ = r;            /* to form our exponent */
  5300.  
  5301.   *p++ = *d++ & 0xff;        /* now do the high order mantissa */
  5302.                 /* strip off the IBM exponent and sign bits */
  5303.   if (mode != SFmode)        /* there are only 2 words in SFmode */
  5304.     {
  5305.       *p++ = *d++;        /* fill in the rest of our mantissa */
  5306.       *p++ = *d++;
  5307.     }
  5308.   *p = *d;
  5309.  
  5310.   if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
  5311.     y[0] = y[E] = 0;
  5312.   else
  5313.     y[E] -= 5 + enormlz (y);    /* now normalise the mantissa */
  5314.                   /* handle change in RADIX */
  5315.   emovo (y, e);
  5316. }
  5317.  
  5318.  
  5319.  
  5320. /* Convert e type to IBM single/double precision.  */
  5321.  
  5322. static void 
  5323. etoibm (x, d, mode)
  5324.      unsigned EMUSHORT *x, *d;
  5325.      enum machine_mode mode;
  5326. {
  5327.   unsigned EMUSHORT xi[NI];
  5328.   EMULONG exp;
  5329.   int rndsav;
  5330.  
  5331.   emovi (x, xi);
  5332.   exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2));    /* adjust exponent for offsets */
  5333.                             /* round off to nearest or even */
  5334.   rndsav = rndprc;
  5335.   rndprc = 56;
  5336.   emdnorm (xi, 0, 0, exp, 64);
  5337.   rndprc = rndsav;
  5338.   toibm (xi, d, mode);
  5339. }
  5340.  
  5341. static void 
  5342. toibm (x, y, mode)
  5343.      unsigned EMUSHORT *x, *y;
  5344.      enum machine_mode mode;
  5345. {
  5346.   unsigned EMUSHORT i;
  5347.   unsigned EMUSHORT *p;
  5348.   int r;
  5349.  
  5350.   p = x;
  5351.   *y = 0;
  5352.   if (*p++)
  5353.     *y = 0x8000;
  5354.   i = *p++;
  5355.   if (i == 0)
  5356.     {
  5357.       *y++ = 0;
  5358.       *y++ = 0;
  5359.       if (mode != SFmode)
  5360.     {
  5361.       *y++ = 0;
  5362.       *y++ = 0;
  5363.     }
  5364.       return;
  5365.     }
  5366.   r = i & 0x3;
  5367.   i >>= 2;
  5368.   if (i > 0x7f)
  5369.     {
  5370.       *y++ |= 0x7fff;
  5371.       *y++ = 0xffff;
  5372.       if (mode != SFmode)
  5373.     {
  5374.       *y++ = 0xffff;
  5375.       *y++ = 0xffff;
  5376.     }
  5377. #ifdef ERANGE
  5378.       errno = ERANGE;
  5379. #endif
  5380.       return;
  5381.     }
  5382.   i &= 0x7f;
  5383.   *y |= (i << 8);
  5384.   eshift (x, r + 5);
  5385.   *y++ |= x[M];
  5386.   *y++ = x[M + 1];
  5387.   if (mode != SFmode)
  5388.     {
  5389.       *y++ = x[M + 2];
  5390.       *y++ = x[M + 3];
  5391.     }
  5392. }
  5393. #endif /* IBM */
  5394.  
  5395. /* Output a binary NaN bit pattern in the target machine's format.  */
  5396.  
  5397. /* If special NaN bit patterns are required, define them in tm.h
  5398.    as arrays of unsigned 16-bit shorts.  Otherwise, use the default
  5399.    patterns here. */
  5400. #ifdef TFMODE_NAN
  5401. TFMODE_NAN;
  5402. #else
  5403. #ifdef MIEEE
  5404. unsigned EMUSHORT TFnan[8] =
  5405.  {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
  5406. #endif
  5407. #ifdef IBMPC
  5408. unsigned EMUSHORT TFnan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
  5409. #endif
  5410. #endif
  5411.  
  5412. #ifdef XFMODE_NAN
  5413. XFMODE_NAN;
  5414. #else
  5415. #ifdef MIEEE
  5416. unsigned EMUSHORT XFnan[6] = {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
  5417. #endif
  5418. #ifdef IBMPC
  5419. unsigned EMUSHORT XFnan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
  5420. #endif
  5421. #endif
  5422.  
  5423. #ifdef DFMODE_NAN
  5424. DFMODE_NAN;
  5425. #else
  5426. #ifdef MIEEE
  5427. unsigned EMUSHORT DFnan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
  5428. #endif
  5429. #ifdef IBMPC
  5430. unsigned EMUSHORT DFnan[4] = {0, 0, 0, 0xfff8};
  5431. #endif
  5432. #endif
  5433.  
  5434. #ifdef SFMODE_NAN
  5435. SFMODE_NAN;
  5436. #else
  5437. #ifdef MIEEE
  5438. unsigned EMUSHORT SFnan[2] = {0x7fff, 0xffff};
  5439. #endif
  5440. #ifdef IBMPC
  5441. unsigned EMUSHORT SFnan[2] = {0, 0xffc0};
  5442. #endif
  5443. #endif
  5444.  
  5445.  
  5446. static void
  5447. make_nan (nan, sign, mode)
  5448.      unsigned EMUSHORT *nan;
  5449.      int sign;
  5450.      enum machine_mode mode;
  5451. {
  5452.   int n;
  5453.   unsigned EMUSHORT *p;
  5454.  
  5455.   switch (mode)
  5456.     {
  5457. /* Possibly the `reserved operand' patterns on a VAX can be
  5458.    used like NaN's, but probably not in the same way as IEEE. */
  5459. #if !defined(DEC) && !defined(IBM)
  5460.     case TFmode:
  5461.       n = 8;
  5462.       p = TFnan;
  5463.       break;
  5464.     case XFmode:
  5465.       n = 6;
  5466.       p = XFnan;
  5467.       break;
  5468.     case DFmode:
  5469.       n = 4;
  5470.       p = DFnan;
  5471.       break;
  5472.     case HFmode:
  5473.     case SFmode:
  5474.       n = 2;
  5475.       p = SFnan;
  5476.       break;
  5477. #endif
  5478.     default:
  5479.       abort ();
  5480.     }
  5481. #ifdef MIEEE
  5482.   *nan++ = (sign << 15) | *p++;
  5483. #endif
  5484.   while (--n != 0)
  5485.     *nan++ = *p++;
  5486. #ifndef MIEEE
  5487.   *nan = (sign << 15) | *p;
  5488. #endif
  5489. }
  5490.  
  5491. /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
  5492.    This is the inverse of the function `etarsingle' invoked by
  5493.    REAL_VALUE_TO_TARGET_SINGLE.  */
  5494.  
  5495. REAL_VALUE_TYPE
  5496. ereal_from_float (f)
  5497.      HOST_WIDE_INT f;
  5498. {
  5499.   REAL_VALUE_TYPE r;
  5500.   unsigned EMUSHORT s[2];
  5501.   unsigned EMUSHORT e[NE];
  5502.  
  5503.   /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
  5504.    This is the inverse operation to what the function `endian' does.  */
  5505. #if FLOAT_WORDS_BIG_ENDIAN
  5506.   s[0] = (unsigned EMUSHORT) (f >> 16);
  5507.   s[1] = (unsigned EMUSHORT) f;
  5508. #else
  5509.   s[0] = (unsigned EMUSHORT) f;
  5510.   s[1] = (unsigned EMUSHORT) (f >> 16);
  5511. #endif
  5512.   /* Convert and promote the target float to E-type. */
  5513.   e24toe (s, e);
  5514.   /* Output E-type to REAL_VALUE_TYPE. */
  5515.   PUT_REAL (e, &r);
  5516.   return r;
  5517. }
  5518.  
  5519.  
  5520. /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
  5521.    This is the inverse of the function `etardouble' invoked by
  5522.    REAL_VALUE_TO_TARGET_DOUBLE.
  5523.  
  5524.    The DFmode is stored as an array of HOST_WIDE_INT in the target's
  5525.    data format, with no holes in the bit packing.  The first element
  5526.    of the input array holds the bits that would come first in the
  5527.    target computer's memory.  */
  5528.  
  5529. REAL_VALUE_TYPE
  5530. ereal_from_double (d)
  5531.      HOST_WIDE_INT d[];
  5532. {
  5533.   REAL_VALUE_TYPE r;
  5534.   unsigned EMUSHORT s[4];
  5535.   unsigned EMUSHORT e[NE];
  5536.  
  5537.   /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
  5538. #if FLOAT_WORDS_BIG_ENDIAN
  5539.   s[0] = (unsigned EMUSHORT) (d[0] >> 16);
  5540.   s[1] = (unsigned EMUSHORT) d[0];
  5541. #if HOST_BITS_PER_WIDE_INT == 32
  5542.   s[2] = (unsigned EMUSHORT) (d[1] >> 16);
  5543.   s[3] = (unsigned EMUSHORT) d[1];
  5544. #else
  5545.   /* In this case the entire target double is contained in the
  5546.      first array element.  The second element of the input is ignored.  */
  5547.   s[2] = (unsigned EMUSHORT) (d[0] >> 48);
  5548.   s[3] = (unsigned EMUSHORT) (d[0] >> 32);
  5549. #endif
  5550. #else
  5551. /* Target float words are little-endian.  */
  5552.   s[0] = (unsigned EMUSHORT) d[0];
  5553.   s[1] = (unsigned EMUSHORT) (d[0] >> 16);
  5554. #if HOST_BITS_PER_WIDE_INT == 32
  5555.   s[2] = (unsigned EMUSHORT) d[1];
  5556.   s[3] = (unsigned EMUSHORT) (d[1] >> 16);
  5557. #else
  5558.   s[2] = (unsigned EMUSHORT) (d[0] >> 32);
  5559.   s[3] = (unsigned EMUSHORT) (d[0] >> 48);
  5560. #endif
  5561. #endif
  5562.   /* Convert target double to E-type. */
  5563.   e53toe (s, e);
  5564.   /* Output E-type to REAL_VALUE_TYPE. */
  5565.   PUT_REAL (e, &r);
  5566.   return r;
  5567. }
  5568.  
  5569.  
  5570. /* Convert target computer unsigned 64-bit integer to e-type.
  5571.    The endian-ness of DImode follows the convention for integers,
  5572.    so we use WORDS_BIG_ENDIAN here, not FLOAT_WORDS_BIG_ENDIAN.  */
  5573.  
  5574. static void
  5575. uditoe (di, e)
  5576.      unsigned EMUSHORT *di;  /* Address of the 64-bit int. */
  5577.      unsigned EMUSHORT *e;
  5578. {
  5579.   unsigned EMUSHORT yi[NI];
  5580.   int k;
  5581.  
  5582.   ecleaz (yi);
  5583. #if WORDS_BIG_ENDIAN
  5584.   for (k = M; k < M + 4; k++)
  5585.     yi[k] = *di++;
  5586. #else
  5587.   for (k = M + 3; k >= M; k--)
  5588.     yi[k] = *di++;
  5589. #endif
  5590.   yi[E] = EXONE + 47;    /* exponent if normalize shift count were 0 */
  5591.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  5592.     ecleaz (yi);        /* it was zero */
  5593.   else
  5594.     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
  5595.   emovo (yi, e);
  5596. }
  5597.  
  5598. /* Convert target computer signed 64-bit integer to e-type. */
  5599.  
  5600. static void
  5601. ditoe (di, e)
  5602.      unsigned EMUSHORT *di;  /* Address of the 64-bit int. */
  5603.      unsigned EMUSHORT *e;
  5604. {
  5605.   unsigned EMULONG acc;
  5606.   unsigned EMUSHORT yi[NI];
  5607.   unsigned EMUSHORT carry;
  5608.   int k, sign;
  5609.  
  5610.   ecleaz (yi);
  5611. #if WORDS_BIG_ENDIAN
  5612.   for (k = M; k < M + 4; k++)
  5613.     yi[k] = *di++;
  5614. #else
  5615.   for (k = M + 3; k >= M; k--)
  5616.     yi[k] = *di++;
  5617. #endif
  5618.   /* Take absolute value */
  5619.   sign = 0;
  5620.   if (yi[M] & 0x8000)
  5621.     {
  5622.       sign = 1;
  5623.       carry = 0;
  5624.       for (k = M + 3; k >= M; k--)
  5625.     {
  5626.       acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
  5627.       yi[k] = acc;
  5628.       carry = 0;
  5629.       if (acc & 0x10000)
  5630.         carry = 1;
  5631.     }
  5632.     }
  5633.   yi[E] = EXONE + 47;    /* exponent if normalize shift count were 0 */
  5634.   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
  5635.     ecleaz (yi);        /* it was zero */
  5636.   else
  5637.     yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
  5638.   emovo (yi, e);
  5639.   if (sign)
  5640.     eneg (e);
  5641. }
  5642.  
  5643.  
  5644. /* Convert e-type to unsigned 64-bit int. */
  5645.  
  5646. static void 
  5647. etoudi (x, i)
  5648.      unsigned EMUSHORT *x;
  5649.      unsigned EMUSHORT *i;
  5650. {
  5651.   unsigned EMUSHORT xi[NI];
  5652.   int j, k;
  5653.  
  5654.   emovi (x, xi);
  5655.   if (xi[0])
  5656.     {
  5657.       xi[M] = 0;
  5658.       goto noshift;
  5659.     }
  5660.   k = (int) xi[E] - (EXONE - 1);
  5661.   if (k <= 0)
  5662.     {
  5663.       for (j = 0; j < 4; j++)
  5664.     *i++ = 0;
  5665.       return;
  5666.     }
  5667.   if (k > 64)
  5668.     {
  5669.       for (j = 0; j < 4; j++)
  5670.     *i++ = 0xffff;
  5671.       if (extra_warnings)
  5672.     warning ("overflow on truncation to integer");
  5673.       return;
  5674.     }
  5675.   if (k > 16)
  5676.     {
  5677.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  5678.      then shift up by 16's.  */
  5679.       j = k - ((k >> 4) << 4);
  5680.       if (j == 0)
  5681.     j = 16;
  5682.       eshift (xi, j);
  5683. #if WORDS_BIG_ENDIAN
  5684.       *i++ = xi[M];
  5685. #else
  5686.       i += 3;
  5687.       *i-- = xi[M];
  5688. #endif
  5689.       k -= j;
  5690.       do
  5691.     {
  5692.       eshup6 (xi);
  5693. #if WORDS_BIG_ENDIAN
  5694.       *i++ = xi[M];
  5695. #else
  5696.       *i-- = xi[M];
  5697. #endif
  5698.     }
  5699.       while ((k -= 16) > 0);
  5700.     }
  5701.   else
  5702.     {
  5703.         /* shift not more than 16 bits */
  5704.       eshift (xi, k);
  5705.  
  5706. noshift:
  5707.  
  5708. #if WORDS_BIG_ENDIAN
  5709.       i += 3;
  5710.       *i-- = xi[M];
  5711.       *i-- = 0;
  5712.       *i-- = 0;
  5713.       *i = 0;
  5714. #else
  5715.       *i++ = xi[M];
  5716.       *i++ = 0;
  5717.       *i++ = 0;
  5718.       *i = 0;
  5719. #endif
  5720.     }
  5721. }
  5722.  
  5723.  
  5724. /* Convert e-type to signed 64-bit int. */
  5725.  
  5726. static void 
  5727. etodi (x, i)
  5728.      unsigned EMUSHORT *x;
  5729.      unsigned EMUSHORT *i;
  5730. {
  5731.   unsigned EMULONG acc;
  5732.   unsigned EMUSHORT xi[NI];
  5733.   unsigned EMUSHORT carry;
  5734.   unsigned EMUSHORT *isave;
  5735.   int j, k;
  5736.  
  5737.   emovi (x, xi);
  5738.   k = (int) xi[E] - (EXONE - 1);
  5739.   if (k <= 0)
  5740.     {
  5741.       for (j = 0; j < 4; j++)
  5742.     *i++ = 0;
  5743.       return;
  5744.     }
  5745.   if (k > 64)
  5746.     {
  5747.       for (j = 0; j < 4; j++)
  5748.     *i++ = 0xffff;
  5749.       if (extra_warnings)
  5750.     warning ("overflow on truncation to integer");
  5751.       return;
  5752.     }
  5753.   isave = i;
  5754.   if (k > 16)
  5755.     {
  5756.       /* Shift more than 16 bits: first shift up k-16 mod 16,
  5757.      then shift up by 16's.  */
  5758.       j = k - ((k >> 4) << 4);
  5759.       if (j == 0)
  5760.     j = 16;
  5761.       eshift (xi, j);
  5762. #if WORDS_BIG_ENDIAN
  5763.       *i++ = xi[M];
  5764. #else
  5765.       i += 3;
  5766.       *i-- = xi[M];
  5767. #endif
  5768.       k -= j;
  5769.       do
  5770.     {
  5771.       eshup6 (xi);
  5772. #if WORDS_BIG_ENDIAN
  5773.       *i++ = xi[M];
  5774. #else
  5775.       *i-- = xi[M];
  5776. #endif
  5777.     }
  5778.       while ((k -= 16) > 0);
  5779.     }
  5780.   else
  5781.     {
  5782.         /* shift not more than 16 bits */
  5783.       eshift (xi, k);
  5784.  
  5785. #if WORDS_BIG_ENDIAN
  5786.       i += 3;
  5787.       *i = xi[M];
  5788.       *i-- = 0;
  5789.       *i-- = 0;
  5790.       *i = 0;
  5791. #else
  5792.       *i++ = xi[M];
  5793.       *i++ = 0;
  5794.       *i++ = 0;
  5795.       *i = 0;
  5796. #endif
  5797.     }
  5798.   /* Negate if negative */
  5799.   if (xi[0])
  5800.     {
  5801.       carry = 0;
  5802. #if WORDS_BIG_ENDIAN
  5803.       isave += 3;
  5804. #endif
  5805.       for (k = 0; k < 4; k++)
  5806.     {
  5807.       acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
  5808. #if WORDS_BIG_ENDIAN
  5809.       *isave-- = acc;
  5810. #else
  5811.       *isave++ = acc;
  5812. #endif
  5813.       carry = 0;
  5814.       if (acc & 0x10000)
  5815.         carry = 1;
  5816.     }
  5817.     }
  5818. }
  5819.  
  5820.  
  5821. /* Longhand square root routine. */
  5822.  
  5823.  
  5824. static int esqinited = 0;
  5825. static unsigned short sqrndbit[NI];
  5826.  
  5827. static void 
  5828. esqrt (x, y)
  5829.      unsigned EMUSHORT *x, *y;
  5830. {
  5831.   unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
  5832.   EMULONG m, exp;
  5833.   int i, j, k, n, nlups;
  5834.  
  5835.   if (esqinited == 0)
  5836.     {
  5837.       ecleaz (sqrndbit);
  5838.       sqrndbit[NI - 2] = 1;
  5839.       esqinited = 1;
  5840.     }
  5841.   /* Check for arg <= 0 */
  5842.   i = ecmp (x, ezero);
  5843.   if (i <= 0)
  5844.     {
  5845.       if (i == -1)
  5846.     {
  5847.       mtherr ("esqrt", DOMAIN);
  5848.       eclear (y);
  5849.     }
  5850.       else
  5851.     emov (x, y);
  5852.       return;
  5853.     }
  5854.  
  5855. #ifdef INFINITY
  5856.   if (eisinf (x))
  5857.     {
  5858.       eclear (y);
  5859.       einfin (y);
  5860.       return;
  5861.     }
  5862. #endif
  5863.   /* Bring in the arg and renormalize if it is denormal. */
  5864.   emovi (x, xx);
  5865.   m = (EMULONG) xx[1];        /* local long word exponent */
  5866.   if (m == 0)
  5867.     m -= enormlz (xx);
  5868.  
  5869.   /* Divide exponent by 2 */
  5870.   m -= 0x3ffe;
  5871.   exp = (unsigned short) ((m / 2) + 0x3ffe);
  5872.  
  5873.   /* Adjust if exponent odd */
  5874.   if ((m & 1) != 0)
  5875.     {
  5876.       if (m > 0)
  5877.     exp += 1;
  5878.       eshdn1 (xx);
  5879.     }
  5880.  
  5881.   ecleaz (sq);
  5882.   ecleaz (num);
  5883.   n = 8;            /* get 8 bits of result per inner loop */
  5884.   nlups = rndprc;
  5885.   j = 0;
  5886.  
  5887.   while (nlups > 0)
  5888.     {
  5889.       /* bring in next word of arg */
  5890.       if (j < NE)
  5891.     num[NI - 1] = xx[j + 3];
  5892.       /* Do additional bit on last outer loop, for roundoff. */
  5893.       if (nlups <= 8)
  5894.     n = nlups + 1;
  5895.       for (i = 0; i < n; i++)
  5896.     {
  5897.       /* Next 2 bits of arg */
  5898.       eshup1 (num);
  5899.       eshup1 (num);
  5900.       /* Shift up answer */
  5901.       eshup1 (sq);
  5902.       /* Make trial divisor */
  5903.       for (k = 0; k < NI; k++)
  5904.         temp[k] = sq[k];
  5905.       eshup1 (temp);
  5906.       eaddm (sqrndbit, temp);
  5907.       /* Subtract and insert answer bit if it goes in */
  5908.       if (ecmpm (temp, num) <= 0)
  5909.         {
  5910.           esubm (temp, num);
  5911.           sq[NI - 2] |= 1;
  5912.         }
  5913.     }
  5914.       nlups -= n;
  5915.       j += 1;
  5916.     }
  5917.  
  5918.   /* Adjust for extra, roundoff loop done. */
  5919.   exp += (NBITS - 1) - rndprc;
  5920.  
  5921.   /* Sticky bit = 1 if the remainder is nonzero. */
  5922.   k = 0;
  5923.   for (i = 3; i < NI; i++)
  5924.     k |= (int) num[i];
  5925.  
  5926.   /* Renormalize and round off. */
  5927.   emdnorm (sq, k, 0, exp, 64);
  5928.   emovo (sq, y);
  5929. }
  5930. #endif /* EMU_NON_COMPILE not defined */
  5931.  
  5932. /* Return the binary precision of the significand for a given
  5933.    floating point mode.  The mode can hold an integer value
  5934.    that many bits wide, without losing any bits.  */
  5935.  
  5936. int
  5937. significand_size (mode)
  5938.      enum machine_mode mode;
  5939. {
  5940.  
  5941. switch (mode)
  5942.   {
  5943.   case SFmode:
  5944.     return 24;
  5945.  
  5946.   case DFmode:
  5947. #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
  5948.     return 53;
  5949. #else
  5950. #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
  5951.     return 56;
  5952. #else
  5953. #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
  5954.     return 56;
  5955. #else
  5956.     abort ();
  5957. #endif
  5958. #endif
  5959. #endif
  5960.  
  5961.   case XFmode:
  5962.     return 64;
  5963.   case TFmode:
  5964.     return 113;
  5965.  
  5966.   default:
  5967.     abort ();
  5968.   }
  5969. }
  5970.